Skip to content

Commit

Permalink
Add MIRS sea-ice product to ioda converter (#1238)
Browse files Browse the repository at this point in the history
#### Description
- Task for adding a new ioda converter for new sea-ice concentration
- For sea-ice, AMSR2 and MIRS product will be used


Partially addressed #1182

---------

Co-authored-by: Guillaume Vernieres <[email protected]>
  • Loading branch information
apchoiCMD and guillaumevernieres authored Jul 30, 2024
1 parent 5583fbb commit 7e09cdc
Show file tree
Hide file tree
Showing 8 changed files with 2,141 additions and 0 deletions.
165 changes: 165 additions & 0 deletions utils/obsproc/IcecMirs2Ioda.h
Original file line number Diff line number Diff line change
@@ -0,0 +1,165 @@
#pragma once

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

#include "eckit/config/LocalConfiguration.h"

#include <Eigen/Dense> // NOLINT

#include "ioda/../../../../core/IodaUtils.h"
#include "ioda/Group.h"
#include "ioda/ObsGroup.h"

#include "oops/util/dateFunctions.h"

#include "NetCDFToIodaConverter.h"

namespace gdasapp {

class IcecMirs2Ioda : public NetCDFToIodaConverter {
public:
explicit IcecMirs2Ioda(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 MIRS" << 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("Scanline").getSize();
int dimySize = ncFile.getDim("Field_of_view").getSize();
int dimQcSize = ncFile.getDim("Qc_dim").getSize();
int nobs = dimxSize * dimySize;
int nqcs = dimxSize * dimySize * dimQcSize;

// 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 non-optional metadata: longitude and latitude
std::vector<float> oneDimLatVal(iodaVars.location_);
ncFile.getVar("Latitude").getVar(oneDimLatVal.data());

std::vector<float> oneDimLonVal(iodaVars.location_);
ncFile.getVar("Longitude").getVar(oneDimLonVal.data());

// Create a vector to hold the full Qc variable
std::vector<int> fullQcVar(nqcs);
ncFile.getVar("Qc").getVar(fullQcVar.data());

// Read Quality Flags as a preQc
// Create a vector to hold the first slice
std::vector<int> oneDimFlagsVal(nobs);
// Extract the first slice
for (int i = 0; i < nobs; ++i) {
oneDimFlagsVal[i] = fullQcVar[i * dimQcSize];
}

// Get Ice_Concentration obs values
std::vector<int16_t> oneDimObsVal(iodaVars.location_);
ncFile.getVar("SIce").getVar(oneDimObsVal.data());

// Read and process the dateTime
std::vector<int64_t> TimeYearVal(dimxSize);
ncFile.getVar("ScanTime_year").getVar(TimeYearVal.data());

std::vector<int64_t> TimeMonthVal(dimxSize);
ncFile.getVar("ScanTime_month").getVar(TimeMonthVal.data());

std::vector<int64_t> TimeDayVal(dimxSize);
ncFile.getVar("ScanTime_dom").getVar(TimeDayVal.data());

std::vector<int64_t> TimeHourVal(dimxSize);
ncFile.getVar("ScanTime_hour").getVar(TimeHourVal.data());

std::vector<int64_t> TimeMinuteVal(dimxSize);
ncFile.getVar("ScanTime_minute").getVar(TimeMinuteVal.data());

std::vector<int64_t> TimeSecondVal(dimxSize);
ncFile.getVar("ScanTime_second").getVar(TimeSecondVal.data());

iodaVars.referenceDate_ = "seconds since 1970-01-01T00:00:00Z";

for (int i = 0; i < dimxSize; i++) {
int year = TimeYearVal[i];
int month = TimeMonthVal[i];
int day = TimeDayVal[i];
int hour = TimeHourVal[i];
int minute = TimeMinuteVal[i];
int second = TimeSecondVal[i];

// Avoid crash util in ioda::convertDtimeToTimeOffsets
if (year == -999 || month == -999 || day == -999 ||
hour == -999 || minute == -999 || second == -999) {
year = month = day = hour = minute = second = 0;
}

// Construct iso8601 string format for each dateTime
std::stringstream ss;
ss << std::setfill('0')
<< std::setw(4) << year << '-'
<< std::setw(2) << month << '-'
<< std::setw(2) << day << 'T'
<< std::setw(2) << hour << ':'
<< std::setw(2) << minute << ':'
<< std::setw(2) << second << 'Z';
std::string formattedDateTime = ss.str();
util::DateTime dateTime(formattedDateTime);

// Set epoch time for MIRS_ICEC
util::DateTime epochDtime("1970-01-01T00:00:00Z");

// Convert Obs DateTime objects to epoch time offsets in seconds
// 0000-00-00T00:00:00Z will be converterd to negative seconds
int64_t timeOffsets
= ioda::convertDtimeToTimeOffsets(epochDtime, {dateTime})[0];

// Update non-optional Eigen arrays
for (int j = 0; j < dimySize; j++) {
int index = i * dimySize + j;
iodaVars.datetime_(index) = timeOffsets;
}
}

// Update non-optional Eigen arrays
for (int i = 0; i < iodaVars.location_; i++) {
iodaVars.longitude_(i) = oneDimLonVal[i];
iodaVars.latitude_(i) = oneDimLatVal[i];
iodaVars.obsVal_(i) = static_cast<float>(oneDimObsVal[i]*0.01);
iodaVars.obsError_(i) = 0.1; // Do something for obs error
iodaVars.preQc_(i) = oneDimFlagsVal[i];
// Store optional metadata, set ocean basins to -999 for now
iodaVars.intMetadata_.row(i) << -999;
}

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

return iodaVars;
};
}; // class IcecMirs2Ioda
} // 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 @@ -8,6 +8,7 @@

#include "../Ghrsst2Ioda.h"
#include "../IcecAmsr2Ioda.h"
#include "../IcecMirs2Ioda.h"
#include "../Rads2Ioda.h"
#include "../RTOFSSalinity.h"
#include "../RTOFSTemperature.h"
Expand Down Expand Up @@ -50,6 +51,9 @@ namespace gdasapp {
} else if (provider == "AMSR2") {
IcecAmsr2Ioda conv2ioda(fullConfig, this->getComm());
conv2ioda.writeToIoda();
} else if (provider == "MIRS") {
IcecMirs2Ioda conv2ioda(fullConfig, this->getComm());
conv2ioda.writeToIoda();
} else if (provider == "VIIRSAOD") {
Viirsaod2Ioda conv2ioda(fullConfig, this->getComm());
conv2ioda.writeToIoda();
Expand Down
9 changes: 9 additions & 0 deletions utils/test/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,7 @@ list( APPEND utils_test_input
testinput/gdas_smap2ioda.yaml
testinput/gdas_smos2ioda.yaml
testinput/gdas_icecamsr2ioda.yaml
testinput/gdas_icecmirs2ioda.yaml
testinput/gdas_viirsaod2ioda.yaml
)

Expand All @@ -19,6 +20,7 @@ set( gdas_utils_test_ref
testref/smap2ioda.test
testref/smos2ioda.test
testref/icecamsr2ioda.test
testref/icecmirs2ioda.test
testref/viirsaod2ioda.test
)

Expand Down Expand Up @@ -145,3 +147,10 @@ ecbuild_add_test( TARGET test_gdasapp_util_icecamsr2ioda
ARGS "../testinput/gdas_icecamsr2ioda.yaml"
LIBS gdas-utils
WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/obsproc)

# Test the MIRS to IODA converter
ecbuild_add_test( TARGET test_gdasapp_util_icecmirs2ioda
COMMAND ${CMAKE_BINARY_DIR}/bin/gdas_obsprovider2ioda.x
ARGS "../testinput/gdas_icecmirs2ioda.yaml"
LIBS gdas-utils
WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/obsproc)
2 changes: 2 additions & 0 deletions utils/test/prepdata.sh
Original file line number Diff line number Diff line change
Expand Up @@ -23,6 +23,8 @@ cdl2nc4 icec_amsr2_north_1.nc4 ${project_source_dir}/testdata/icec_amsr2_north_1
cdl2nc4 icec_amsr2_north_2.nc4 ${project_source_dir}/testdata/icec_amsr2_north_2.cdl
cdl2nc4 icec_amsr2_south_1.nc4 ${project_source_dir}/testdata/icec_amsr2_south_1.cdl
cdl2nc4 icec_amsr2_south_2.nc4 ${project_source_dir}/testdata/icec_amsr2_south_2.cdl
cdl2nc4 icec_mirs_snpp_1.nc4 ${project_source_dir}/testdata/icec_mirs_snpp_1.cdl
cdl2nc4 icec_mirs_snpp_2.nc4 ${project_source_dir}/testdata/icec_mirs_snpp_2.cdl
cdl2nc4 sss_smap_1.nc4 ${project_source_dir}/testdata/sss_smap_1.cdl
cdl2nc4 sss_smap_2.nc4 ${project_source_dir}/testdata/sss_smap_2.cdl
cdl2nc4 sss_smos_1.nc4 ${project_source_dir}/testdata/sss_smos_1.cdl
Expand Down
Loading

0 comments on commit 7e09cdc

Please sign in to comment.