Merge branch 'develop' into feature/aero-jcb

modulefiles/GDAS/hercules.gnu.lua

@@ -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")

test/marine/testinput/bufr2ioda_insitu_profile_tropical_2019010700.yaml.in

@@ -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
 

utils/obsproc/IcecAbi2Ioda.h

@@ -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

utils/obsproc/applications/gdas_obsprovider2ioda.h

@@ -7,6 +7,7 @@
 #include "oops/runs/Application.h"
 
 #include "../Ghrsst2Ioda.h"
+#include "../IcecAbi2Ioda.h"
 #include "../IcecAmsr2Ioda.h"
 #include "../IcecJpssrr2Ioda.h"
 #include "../IcecMirs2Ioda.h"
@@ -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();

utils/obsproc/util.h

@@ -1,6 +1,7 @@
 #pragma once
 
 #include <iostream>
+#include <limits>
 #include <map>
 #include <netcdf> // NOLINT (using C API)
 #include <string>
@@ -254,5 +255,69 @@ namespace gdasapp {
  }
  };
  } // namespace iodavars
+
+ // TODO(Mindo): To move below as a private method to the iceabi2ioda class
+ namespace utils {
+
+ // Calculate latitude and longitude from GOES ABI fixed grid projection data
+ // GOES ABI fixed grid projection is a map projection relative to the GOES satellite
+ // Units: latitude in °N (°S < 0), longitude in °E (°W < 0)
+ // See GOES-R Product User Guide (PUG) Volume 5 (L2 products) Section 4.2.8 (p58)
+ void abiToGeoLoc(
+ const std::vector<std::vector<double>>& x_coordinate_2d,
+ const std::vector<std::vector<double>>& y_coordinate_2d,
+ double lon_origin,
+ double H,
+ double r_eq,
+ double r_pol,
+ std::vector<std::vector<double>>& abi_lat,
+ std::vector<std::vector<double>>& abi_lon
+ ) {
+ int sizeX = x_coordinate_2d[0].size();
+ int sizeY = x_coordinate_2d.size();
+
+ double lambda_0 = (lon_origin * M_PI) / 180.0;
+
+ abi_lat.resize(sizeY, std::vector<double>(sizeX));
+ abi_lon.resize(sizeY, std::vector<double>(sizeX));
+
+ for (int i = 0; i < sizeY; ++i) {
+ for (int j = 0; j < sizeX; ++j) {
+ double x = x_coordinate_2d[i][j];
+ double y = y_coordinate_2d[i][j];
+
+ // Cache sin(x), cos(x), sin(y), and cos(y)
+ double sin_x = std::sin(x);
+ double cos_x = std::cos(x);
+ double sin_y = std::sin(y);
+ double cos_y = std::cos(y);
+
+ double a_var = std::pow(sin_x, 2.0)
+ + std::pow(cos_x, 2.0) * (std::pow(cos_y, 2.0)
+ + ((r_eq * r_eq) / (r_pol * r_pol)) * std::pow(sin_y, 2.0));
+ double b_var = -2.0 * H * cos_x * cos_y;
+ double c_var = (H * H) - (r_eq * r_eq);
+ double discriminant = (b_var * b_var) - (4.0 * a_var * c_var);
+
+ // Check if discriminant is strictly positive
+ if (discriminant > 0) {
+ double r_s = (-b_var - std::sqrt(discriminant)) / (2.0 * a_var);
+ double s_x = r_s * cos_x * cos_y;
+ double s_y = -r_s * sin_x;
+ double s_z = r_s * cos_x * sin_y;
+
+ abi_lat[i][j] = (180.0 / M_PI) * (std::atan(((r_eq * r_eq) / (r_pol * r_pol))
+ * (s_z / std::sqrt(((H - s_x) * (H - s_x)) + (s_y * s_y)))));
+ abi_lon[i][j] = (lambda_0 - std::atan(s_y / (H - s_x))) * (180.0 / M_PI);
+ } else {
+ // Handle invalid values
+ // Set latitude and longitude to NaN if discriminant <= 0)
+ abi_lat[i][j] = std::numeric_limits<double>::quiet_NaN();
+ abi_lon[i][j] = std::numeric_limits<double>::quiet_NaN();
+ }
+ }
+ }
+ } // void
+ } // namespace utils
  } // namespace obsproc
 }; // namespace gdasapp