Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add MIRS sea-ice product to ioda converter #1238

Merged
merged 3 commits into from
Jul 30, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
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);
guillaumevernieres marked this conversation as resolved.
Show resolved Hide resolved

// 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
Loading