Skip to content

Commit

Permalink
Merge branch 'develop' at f1222ec into patch/gwci
Browse files Browse the repository at this point in the history
  • Loading branch information
RussTreadon-NOAA committed Oct 23, 2024
2 parents 0325836 + f1222ec commit 764f58c
Show file tree
Hide file tree
Showing 14 changed files with 1,723 additions and 5 deletions.
93 changes: 93 additions & 0 deletions modulefiles/GDAS/hercules.gnu.lua
Original file line number Diff line number Diff line change
@@ -0,0 +1,93 @@
help([[
Load environment for running the GDAS application with gnu compilers and MPI.
]])

local pkgName = myModuleName()
local pkgVersion = myModuleVersion()
local pkgNameVer = myModuleFullName()

prepend_path("MODULEPATH", '/work/noaa/epic/role-epic/spack-stack/hercules/modulefiles')
prepend_path("MODULEPATH", '/work/noaa/epic/role-epic/spack-stack/hercules/spack-stack-1.7.0/envs/ue-gcc/install/modulefiles/Core')
prepend_path("MODULEPATH", '/work2/noaa/da/python/opt/modulefiles/stack')


---- below two lines get us access to the spack-stack modules
load("stack-gcc/12.2.0")
load("stack-openmpi/4.1.6")

load("cmake/3.23.1")
load("curl/8.4.0")
load("zlib/1.2.13")
load("git/2.31.1")
--load("pkg-config/0.27.1")
load("hdf5/1.14.3")
load("parallel-netcdf/1.12.3")
load("netcdf-c/4.9.2")
load("nccmp/1.9.0.1")
load("netcdf-fortran/4.6.1")
load("nco/5.1.6")
load("parallelio/2.6.2")
load("wget/1.21.1")
load("boost/1.84.0")
load("bufr/12.0.1")
load("git-lfs/3.1.2")
load("ecbuild/3.7.2")
load("openjpeg/2.3.1")
load("eccodes/2.33.0")
load("eigen/3.4.0")
load("openblas/0.3.24")
load("eckit/1.24.5")
load("fftw/3.3.10")
load("fckit/0.11.0")
load("fiat/1.2.0")
load("ectrans/1.2.0")
load("fms/2023.04")
load("esmf/8.6.1")
load("atlas/0.36.0")
load("sp/2.5.0")
load("gsl-lite/0.37.0")
load("libjpeg/2.1.0")
load("krb5/1.20.1")
load("libtirpc/1.3.3")
load("hdf/4.2.15")
load("jedi-cmake/1.4.0")
load("libpng/1.6.37")
load("libxt/1.3.0")
load("libxmu/1.1.4")
load("libxpm/3.5.17")
load("libxaw/1.0.15")
load("udunits/2.2.28")
load("ncview/2.1.9")
load("netcdf-cxx4/4.3.1")
load("py-pybind11/2.11.0")
--load("crtm/v2.4_jedi")
load("contrib/0.1")
load("noaatools/3.1")
load("rocoto/1.3.7")

load("hpc/1.2.0")
unload("python/3.10.13")
unload("py-numpy/1.22.3")
load("miniconda3/4.6.14")
load("gdasapp/1.0.0")
-- below is a hack because of cmake finding the wrong python...
setenv("CONDA_PREFIX", "/work2/noaa/da/python/opt/core/miniconda3/4.6.14/envs/gdasapp/")

setenv("CC","mpicc")
setenv("FC","mpifort")
setenv("CXX","mpicxx")
local mpiexec = '/opt/slurm/bin/srun'
local mpinproc = '-n'
setenv('MPIEXEC_EXEC', mpiexec)
setenv('MPIEXEC_NPROC', mpinproc)

setenv("CRTM_FIX","/work2/noaa/da/role-da/GDASApp/fix/crtm/2.4.0")
setenv("GDASAPP_TESTDATA","/work2/noaa/da/role-da/GDASApp/testdata")
setenv("GDASAPP_UNIT_TEST_DATA_PATH", "/work2/noaa/da/role-da/GDASApp/unittestdata")

execute{cmd="ulimit -s unlimited",modeA={"load"}}

whatis("Name: ".. pkgName)
whatis("Version: ".. pkgVersion)
whatis("Category: GDASApp")
whatis("Description: Load all libraries needed for GDASApp")
2 changes: 1 addition & 1 deletion parm/aero/obs/config/viirs_n20_aod.yaml.j2
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@
obsdatain:
engine:
type: H5File
obsfile: "{{ DATA }}/obs/{{ OPREFIX }}viirs_n20.{{ current_cycle | to_YMDH }}.nc4"
obsfile: "{{ DATA }}/obs/{{ OPREFIX }}viirs_n20_aod.tm00.nc"
obsdataout:
engine:
type: H5File
Expand Down
2 changes: 1 addition & 1 deletion parm/aero/obs/config/viirs_n21_aod.yaml.j2
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@
obsdatain:
engine:
type: H5File
obsfile: "{{ DATA }}/obs/{{ OPREFIX }}viirs_n21.{{ current_cycle | to_YMDH }}.nc4"
obsfile: "{{ DATA }}/obs/{{ OPREFIX }}viirs_n21_aod.tm00.nc"
obsdataout:
engine:
type: H5File
Expand Down
2 changes: 1 addition & 1 deletion parm/aero/obs/config/viirs_npp_aod.yaml.j2
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@
obsdatain:
engine:
type: H5File
obsfile: "{{ DATA }}/obs/{{ OPREFIX }}viirs_npp.{{ current_cycle | to_YMDH }}.nc4"
obsfile: "{{ DATA }}/obs/{{ OPREFIX }}viirs_npp_aod.tm00.nc"
obsdataout:
engine:
type: H5File
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -2,12 +2,12 @@
data_format: dbuoy
subsets: dbuoy
source: NCEP data tank
data_type: drifter
data_type: tropical
cycle_type: gdas
cycle_datetime: '2019010700'
dump_directory: __BUFRINPUTDIR__
ioda_directory: __IODAOUTPUTDIR__
ocean_basin: __OCEANBASIN__
data_description: 6-hrly in situ drifter profiles
data_description: 6-hrly in situ tropical mooring profiles
data_provider: U.S. NOAA

177 changes: 177 additions & 0 deletions utils/obsproc/IcecAbi2Ioda.h
Original file line number Diff line number Diff line change
@@ -0,0 +1,177 @@
#pragma once

#include <ctime>
#include <iomanip>
#include <iostream>
#include <map>
#include <netcdf> // NOLINT (using C API)
#include <string>
#include <vector>

#include "eckit/config/LocalConfiguration.h"

#include <Eigen/Dense> // NOLINT

#include "ioda/../../../../core/IodaUtils.h" // TODO(All): Use a better way in all converters
#include "ioda/Group.h"
#include "ioda/ObsGroup.h"

#include "oops/util/dateFunctions.h"

#include "NetCDFToIodaConverter.h"

namespace gdasapp {

class IcecAbi2Ioda : public NetCDFToIodaConverter {
public:
explicit IcecAbi2Ioda(const eckit::Configuration & fullConfig, const eckit::mpi::Comm & comm)
: NetCDFToIodaConverter(fullConfig, comm) {
variable_ = "seaIceFraction";
}

// Read netcdf file and populate iodaVars
gdasapp::obsproc::iodavars::IodaVars providerToIodaVars(const std::string fileName) final {
oops::Log::info() << "Processing files provided by the ABI" << std::endl;

// Open the NetCDF file in read-only mode
netCDF::NcFile ncFile(fileName, netCDF::NcFile::read);
oops::Log::info() << "Reading... " << fileName << std::endl;

// Get the number of obs in the file
int dimxSize = ncFile.getDim("x").getSize();
int dimySize = ncFile.getDim("y").getSize();
int nobs = dimxSize * dimySize;

// Set the int metadata names
std::vector<std::string> intMetadataNames = {"oceanBasin"};

// Set the float metadata name
std::vector<std::string> floatMetadataNames = {};

// Create instance of iodaVars object
gdasapp::obsproc::iodavars::IodaVars iodaVars(nobs, floatMetadataNames, intMetadataNames);

oops::Log::debug() << "--- iodaVars.location_: " << iodaVars.location_ << std::endl;

// Read in GOES ABI fixed grid projection variables and constants
std::vector<double> y_coordinate_1d(dimySize);
ncFile.getVar("y").getVar(y_coordinate_1d.data());
float yOffSet;
ncFile.getVar("y").getAtt("add_offset").getValues(&yOffSet);
float yScaleFactor;
ncFile.getVar("y").getAtt("scale_factor").getValues(&yScaleFactor);
// Apply the scale factor and add offset to the raw data
for (auto& yval : y_coordinate_1d) {
yval = yval * yScaleFactor + yOffSet; // N-S elevation angle in radians
}

std::vector<double> x_coordinate_1d(dimxSize);
ncFile.getVar("x").getVar(x_coordinate_1d.data());
float xOffSet;
ncFile.getVar("x").getAtt("add_offset").getValues(&xOffSet);
float xScaleFactor;
ncFile.getVar("x").getAtt("scale_factor").getValues(&xScaleFactor);
// Apply the scale factor and add offset to the raw data
for (auto& xval : x_coordinate_1d) {
xval = xval * xScaleFactor + xOffSet; // E-W scanning angle in radians
}

// Create 2D arrays (meshgrid equivalent)
std::vector<std::vector<double>> x_coordinate_2d(dimySize, std::vector<double>(dimxSize));
std::vector<std::vector<double>> y_coordinate_2d(dimySize, std::vector<double>(dimxSize));
std::vector<std::vector<double>> abi_lon;
std::vector<std::vector<double>> abi_lat;

// Create 2D coordinate matrices from 1D coordinate vectors
for (int i = 0; i < dimySize; ++i) {
for (int j = 0; j < dimxSize; ++j) {
x_coordinate_2d[i][j] = x_coordinate_1d[j];
y_coordinate_2d[i][j] = y_coordinate_1d[i];
}
}

// Retrieve the attributes
double lon_origin;
ncFile.getVar("goes_imager_projection").getAtt("longitude_of_projection_origin")
.getValues(&lon_origin);
double perspective_point_height;
ncFile.getVar("goes_imager_projection").getAtt("perspective_point_height")
.getValues(&perspective_point_height);
double r_eq;
ncFile.getVar("goes_imager_projection").getAtt("semi_major_axis").getValues(&r_eq);
double r_pol;
ncFile.getVar("goes_imager_projection").getAtt("semi_minor_axis").getValues(&r_pol);

// Calculate H = Satellite height from center of earth(m)
double H = perspective_point_height + r_eq;

// Calculate Latitude and Longitude from GOES Imager Projection
// for details of calculations in util.h
gdasapp::obsproc::utils::abiToGeoLoc(
x_coordinate_2d,
y_coordinate_2d,
lon_origin,
H,
r_eq,
r_pol,
abi_lat,
abi_lon);

// Store real number of lat and lon into eigen arrays
int loc(0);
for (int i = 0; i < dimySize; i++) {
for (int j = 0; j < dimxSize; j++) {
iodaVars.longitude_(loc) = std::real(abi_lon[i][j]);
iodaVars.latitude_(loc) = std::real(abi_lat[i][j]);
loc += 1;
}
}

// Read Quality Flags as a preQc
std::vector<uint16_t> fullQcFlagsVar(iodaVars.location_);
ncFile.getVar("DQF").getVar(fullQcFlagsVar.data());

// Get Ice_Concentration obs values
std::vector<uint16_t> IcecObsVal(iodaVars.location_);
ncFile.getVar("IceConc").getVar(IcecObsVal.data());
float IcecOffSet;
ncFile.getVar("IceConc").getAtt("add_offset").getValues(&IcecOffSet);
float IcecScaleFactor;
ncFile.getVar("IceConc").getAtt("scale_factor").getValues(&IcecScaleFactor);

// TODO(All): Think how we will be acle to use Temp later
// Get Ice_Temp obs values
std::vector<uint16_t> IcecTempObsVal(iodaVars.location_);
ncFile.getVar("Temp").getVar(IcecTempObsVal.data()); // Kelvin
float IcecTempOffSet;
ncFile.getVar("Temp").getAtt("add_offset").getValues(&IcecTempOffSet);
float IcecTempScaleFactor;
ncFile.getVar("Temp").getAtt("scale_factor").getValues(&IcecTempScaleFactor);

// Read the dateTime
double TimeVal;
ncFile.getVar("t").getVar(&TimeVal);

iodaVars.referenceDate_ = "seconds since 2000-01-01T12:00:00Z"; // 12Z

// Update Eigen arrays
for (int i = 0; i < iodaVars.location_; i++) {
iodaVars.obsVal_(i)
= static_cast<float>((IcecObsVal[i] * IcecScaleFactor + IcecOffSet)*0.01);
iodaVars.obsError_(i) = 0.1; // Do something for obs error
iodaVars.preQc_(i) = fullQcFlagsVar[i];
// Store optional metadata, set ocean basins to -999 for now
iodaVars.intMetadata_.row(i) << -999;
iodaVars.datetime_(i) = TimeVal;
}

// basic test for iodaVars.trim
Eigen::Array<bool, Eigen::Dynamic, 1> mask =
((iodaVars.obsVal_ >= 0.0 && iodaVars.obsVal_ <= 1.0) &&
(iodaVars.latitude_ <= -40.0 || iodaVars.latitude_ >= 40.0));
iodaVars.trim(mask);

return iodaVars;
};
}; // class IcecAbi2Ioda
} // namespace gdasapp
4 changes: 4 additions & 0 deletions utils/obsproc/applications/gdas_obsprovider2ioda.h
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,7 @@
#include "oops/runs/Application.h"

#include "../Ghrsst2Ioda.h"
#include "../IcecAbi2Ioda.h"
#include "../IcecAmsr2Ioda.h"
#include "../IcecJpssrr2Ioda.h"
#include "../IcecMirs2Ioda.h"
Expand Down Expand Up @@ -49,6 +50,9 @@ namespace gdasapp {
} else if (provider == "SMOS") {
Smos2Ioda conv2ioda(fullConfig, this->getComm());
conv2ioda.writeToIoda();
} else if (provider == "ABI") {
IcecAbi2Ioda conv2ioda(fullConfig, this->getComm());
conv2ioda.writeToIoda();
} else if (provider == "AMSR2") {
IcecAmsr2Ioda conv2ioda(fullConfig, this->getComm());
conv2ioda.writeToIoda();
Expand Down
Loading

0 comments on commit 764f58c

Please sign in to comment.