diff --git a/src/Calibrator/Altitude.cpp b/src/Calibrator/Altitude.cpp index 99617e8d..c6351f27 100644 --- a/src/Calibrator/Altitude.cpp +++ b/src/Calibrator/Altitude.cpp @@ -30,6 +30,6 @@ bool CalibratorAltitude::calibrateCore(File& iFile, const ParameterFile* iParame std::string CalibratorAltitude::description() { std::stringstream ss; - ss << Util::formatDescription("-c altitude", "Changes the altitudes in the file to the altitudes in the parameter file.") << std::endl; + ss << Util::formatDescription("-c altitude", "Changes the altitudes in the file to the altitudes in the parameter file. If the file does not have an altitude field, then a new one is not created.") << std::endl; return ss.str(); } diff --git a/src/File/Ec.cpp b/src/File/Ec.cpp index 5352365d..c8cc2a5b 100644 --- a/src/File/Ec.cpp +++ b/src/File/Ec.cpp @@ -151,6 +151,7 @@ void FileEc::writeCore(std::vector iVariables) { writeTimes(); writeReferenceTime(); + writeAltitude(); for(int v = 0; v < iVariables.size(); v++) { Variable::Type varType = iVariables[v]; std::string variable = getVariableName(varType); @@ -390,3 +391,57 @@ std::string FileEc::description() { ss << Util::formatDescription("type=ec", "ECMWF ensemble file") << std::endl; return ss.str(); } + +void FileEc::writeAltitude() const { + if(!hasVar("altitude")) { + return; + } + int vElev = getVar("altitude"); + int numDims = getNumDims(vElev); + if(numDims == 1) { + Util::error("Cannot write altitude when the variable only has one dimension"); + } + + int N = getNumDims(vElev); + long count[N]; + int size = 1; + int indexLat = Util::MV; + int indexLon = Util::MV; + int dims[N]; + nc_inq_vardimid(mFile, vElev, dims); + for(int i = 0; i < N; i++) { + if(dims[i] == getLatDim()) { + count[i] = getNumLat(); + size *= count[i]; + indexLat = i; + } + else if(dims[i] == getLonDim()) { + count[i] = getNumLon(); + size *= count[i]; + indexLon = i; + } + else { + count[i] = 1; + } + } + float MV = getMissingValue(vElev); + float* values = new float[size]; + vec2 elevs = getElevs(); + for(int i = 0; i < getNumLat(); i++) { + for(int j = 0; j < getNumLon(); j++) { + float elev = elevs[i][j]; + if(Util::isValid(MV) && !Util::isValid(elev)) + elev = MV; + // Latitude dimension is ordered first + if(indexLat < indexLon) { + values[i*getNumLon() + j] = elev; + } + // Longitude dimension is ordered first + else { + values[j*getNumLat() + i] = elev; + } + } + } + nc_put_var_float(mFile, vElev, values); + delete[] values; +} diff --git a/src/File/Ec.h b/src/File/Ec.h index a646bb34..0cfe6c70 100644 --- a/src/File/Ec.h +++ b/src/File/Ec.h @@ -29,6 +29,7 @@ class FileEc : public FileNetcdf { std::vector mTimes; vec2 getGridValues(int iVariable) const; + void writeAltitude() const; int getLatDim() const; int getLonDim() const; int getLatVar() const; diff --git a/src/Testing/CalibratorAltitude.cpp b/src/Testing/CalibratorAltitude.cpp index 63250749..8f4bf26e 100644 --- a/src/Testing/CalibratorAltitude.cpp +++ b/src/Testing/CalibratorAltitude.cpp @@ -16,7 +16,7 @@ namespace { virtual void TearDown() { } }; - TEST_F(TestCalibratorAltitude, 10x10) { + TEST_F(TestCalibratorAltitude, arome) { FileArome from("testing/files/10x10.nc"); ParameterFileNetcdf par(Options("file=testing/files/10x10_param_zero_altitude.nc")); CalibratorAltitude cal = CalibratorAltitude(Options()); @@ -36,6 +36,26 @@ namespace { EXPECT_FLOAT_EQ(304, (*after)(5,9,0)); EXPECT_FLOAT_EQ(320, (*after)(0,9,0)); } + TEST_F(TestCalibratorAltitude, ec) { + FileEc from("testing/files/10x10_ec.nc"); + ParameterFileNetcdf par(Options("file=testing/files/10x10_param_zero_altitude.nc")); + CalibratorAltitude cal = CalibratorAltitude(Options()); + + cal.calibrate(from, &par); + + // Elevations should all be 0 + vec2 elevs = from.getElevs(); + for(int i = 0; i < from.getNumLat(); i++) { + for(int j = 0; j < from.getNumLon(); j++) { + EXPECT_FLOAT_EQ(0, elevs[i][j]); + } + } + // Shouldn't have changed anything else + FieldPtr after = from.getField(Variable::T, 0); + EXPECT_FLOAT_EQ(301, (*after)(5,2,0)); + EXPECT_FLOAT_EQ(304, (*after)(5,9,0)); + EXPECT_FLOAT_EQ(320, (*after)(0,9,0)); + } // Should not calibrate when a parameter file with no locaions is used TEST_F(TestCalibratorAltitude, locationIndependent) { ::testing::FLAGS_gtest_death_test_style = "threadsafe"; diff --git a/testing/files/10x10_ec.nc b/testing/files/10x10_ec.nc new file mode 100644 index 00000000..6f39662c Binary files /dev/null and b/testing/files/10x10_ec.nc differ