Merge branch 'mask'

debian/changelog

@@ -1,3 +1,9 @@
+gridpp (0.2.7) precise; urgency=low
+
+  * Adds -c mask
+
+ -- Thomas Nipen <thomas.nipen@met.no>  Thu, 17 Nov 2016 15:36:47 +0100
+
 gridpp (0.2.6) precise; urgency=low
 
   * Adds precipitation rate variable

src/Calibrator/Calibrator.cpp

@@ -157,6 +157,15 @@ Calibrator* Calibrator::getScheme(std::string iName, const Options& iOptions) {
 
       return c;
    }
+   else if(iName == "mask") {
+      std::string variable;
+      if(!iOptions.getValue("variable", variable)) {
+         Util::error("Calibrator 'mask' needs variable");
+      }
+      CalibratorMask* c = new CalibratorMask(Variable::getType(variable), iOptions);
+
+      return c;
+   }
    else {
       Util::error("Could not instantiate calibrator with name '" + iName + "'");
       return NULL;
@@ -213,6 +222,7 @@ std::string Calibrator::getDescriptions() {
    ss << CalibratorGaussian::description() << std::endl;
    ss << CalibratorKriging::description() << std::endl;
    ss << CalibratorNeighbourhood::description() << std::endl;
+   ss << CalibratorMask::description() << std::endl;
    ss << CalibratorPhase::description() << std::endl;
    ss << CalibratorQc::description() << std::endl;
    ss << CalibratorQnh::description() << std::endl;

src/Calibrator/Calibrator.h

@@ -52,6 +52,7 @@ class Calibrator : public Scheme {
 #include "Gaussian.h"
 #include "Kriging.h"
 #include "Neighbourhood.h"
+#include "Mask.h"
 #include "Phase.h"
 #include "Qc.h"
 #include "Qnh.h"

src/Calibrator/Mask.cpp

@@ -0,0 +1,82 @@
+#include "Mask.h"
+#include <cmath>
+#include "../Util.h"
+#include "../File/File.h"
+#include "../ParameterFile/ParameterFile.h"
+CalibratorMask::CalibratorMask(Variable::Type iVariable, const Options& iOptions) :
+      Calibrator(iOptions),
+      mVariable(iVariable),
+      mUseNearestOnly(false),
+      mKeep(true) {
+   iOptions.getValue("keep", mKeep);
+}
+bool CalibratorMask::calibrateCore(File& iFile, const ParameterFile* iParameterFile) const {
+   if(!iParameterFile->isLocationDependent()) {
+      Util::error("Cannot use a location independent parameter file in mask calibrator");
+   }
+   if(!iParameterFile->isFixedSize()) {
+      Util::error("Cannot use a parameter file without a constant number of parameters");
+   }
+   int nLat = iFile.getNumLat();
+   int nLon = iFile.getNumLon();
+   int nEns = iFile.getNumEns();
+   int nTime = iFile.getNumTime();
+   vec2 lats = iFile.getLats();
+   vec2 lons = iFile.getLons();
+   vec2 elevs = iFile.getElevs();
+   vec2 isWithinRadius;
+
+   // Stores if a certain gridpoint is within some radius of influence
+   isWithinRadius.resize(nLat);
+   for(int i = 0; i < nLat; i++) {
+      isWithinRadius[i].resize(nLon, 0);
+   }
+
+   for(int t = 0; t < nTime; t++) {
+      FieldPtr field = iFile.getField(mVariable, t);
+      #pragma omp parallel for
+      for(int i = 0; i < nLat; i++) {
+         for(int j = 0; j < nLon; j++) {
+            if(t == 0 || iParameterFile->isTimeDependent()) {
+               Location loc(Util::MV, Util::MV, Util::MV);
+               Location currLocation(lats[i][j], lons[i][j], elevs[i][j]);
+               if(mUseNearestOnly) {
+                  iParameterFile->getNearestLocation(t, currLocation, loc);
+                  Parameters parameters = iParameterFile->getParameters(t, Location(lats[i][j], lons[i][j], elevs[i][j]));
+                  float dist = currLocation.getDistance(loc);
+                  if(Util::isValid(dist) && dist < parameters[0]) {
+                     isWithinRadius[i][j] = 1;
+                  }
+               }
+               else {
+                  std::vector<Location> locations = iParameterFile->getLocations();
+                  for(int k = 0; k < locations.size(); k++) {
+                     Parameters parameters = iParameterFile->getParameters(t, locations[k]);
+                     float dist = currLocation.getDistance(locations[k]);
+                     if(Util::isValid(dist) && dist < parameters[0]) {
+                        isWithinRadius[i][j] = 1;
+                     }
+                  }
+               }
+            }
+            // Remove the point if keep and we are not within the radius or we don't keep and are
+            // inside
+            bool remove = mKeep != isWithinRadius[i][j];
+            if(remove) {
+               for(int e = 0; e < nEns; e++) {
+                  (*field)(i,j,e) = Util::MV;
+               }
+            }
+         }
+      }
+   }
+   iFile.setElevs(elevs);
+   return true;
+}
+
+std::string CalibratorMask::description() {
+   std::stringstream ss;
+   ss << Util::formatDescription("-c mask", "Allows removal of gridpoints that are not within (or alternatively within) a radius of a set of points provided in a parameter file. The removed gridpoints are set to missing. The parameter file must be spatial and conain one parameter representing the radius in meters. Each gridpoint is checked against each parameter point to see if the gridpoint is within its radius.") << std::endl;
+   ss << Util::formatDescription("   keep=1", "If 1, then all gridpoints within the radius of any parameter point will be kept. If 0, then all gridpoints within the radius of any parameter point will be removed.") << std::endl;
+   return ss.str();
+}

src/Calibrator/Mask.h

@@ -0,0 +1,19 @@
+#ifndef CALIBRATOR_MASK_H
+#define CALIBRATOR_MASK_H
+#include "Calibrator.h"
+class ParameterFile;
+class Parameters;
+
+//! Sets gridpoints that are far away from parameter points to missing
+class CalibratorMask : public Calibrator {
+   public:
+      CalibratorMask(Variable::Type iVariable, const Options& iOptions);
+      static std::string description();
+      std::string name() const {return "mask";};
+   private:
+      bool calibrateCore(File& iFile, const ParameterFile* iParameterFile) const;
+      Variable::Type mVariable;
+      bool mUseNearestOnly;
+      bool mKeep;
+};
+#endif

src/File/Arome.cpp

@@ -534,9 +534,9 @@ bool FileArome::isValid(std::string iFilename) {
 std::string FileArome::description() {
    std::stringstream ss;
    ss << Util::formatDescription("type=arome", "AROME file") << std::endl;
-   ss << Util::formatDescription("lat=latitude", "Name of the variable representing latitudes") << std::endl;
-   ss << Util::formatDescription("lon=longitude", "Name of the variable representing longitudes") << std::endl;
-   ss << Util::formatDescription("x=undef", "Name of dimension in the x-direction. If unspecified, the name is auto-detected.") << std::endl;
-   ss << Util::formatDescription("y=undef", "Name of dimension in the y-direction. If unspecified, the name is auto-detected.") << std::endl;
+   ss << Util::formatDescription("   lat=latitude", "Name of the variable representing latitudes") << std::endl;
+   ss << Util::formatDescription("   lon=longitude", "Name of the variable representing longitudes") << std::endl;
+   ss << Util::formatDescription("   x=undef", "Name of dimension in the x-direction. If unspecified, the name is auto-detected.") << std::endl;
+   ss << Util::formatDescription("   y=undef", "Name of dimension in the y-direction. If unspecified, the name is auto-detected.") << std::endl;
    return ss.str();
 }

src/Testing/CalibratorMask.cpp

@@ -0,0 +1,78 @@
+#include "../File/Fake.h"
+#include "../Util.h"
+#include "../ParameterFile/ParameterFile.h"
+#include "../Calibrator/Mask.h"
+#include <gtest/gtest.h>
+
+namespace {
+   class TestCalibratorMask : public ::testing::Test {
+      protected:
+         TestCalibratorMask() {
+         }
+         virtual ~TestCalibratorMask() {
+         }
+         virtual void SetUp() {
+         }
+         virtual void TearDown() {
+         }
+         std::vector<float> getVector(float iArray[]) {
+            return std::vector<float>(iArray, iArray + sizeof(iArray)/sizeof(float));
+         }
+   };
+   // Mask out values
+   TEST_F(TestCalibratorMask, 10x10_mask_out) {
+      FileArome from("testing/files/10x10.nc");
+      // One point at 3,5 with 223 km radius and one point at 4,6 with 336 km radius
+      ParameterFileText par(Options("file=testing/files/mask0.txt spatial=1"));
+      CalibratorMask cal = CalibratorMask(Variable::T, Options("keep=0"));
+
+      cal.calibrate(from, &par);
+      FieldPtr after = from.getField(Variable::T, 0);
+      ASSERT_EQ(10, after->getNumLat());
+      ASSERT_EQ(10, after->getNumLon());
+      ASSERT_EQ(1,  after->getNumEns());
+
+      EXPECT_FLOAT_EQ(301, (*after)(5,2,0));
+      EXPECT_FLOAT_EQ(Util::MV, (*after)(3,5,0));
+      EXPECT_FLOAT_EQ(Util::MV, (*after)(3,3,0));
+      EXPECT_FLOAT_EQ(Util::MV, (*after)(2,5,0));
+      EXPECT_FLOAT_EQ(Util::MV, (*after)(4,9,0));
+      EXPECT_FLOAT_EQ(302, (*after)(2,3,0));
+      EXPECT_FLOAT_EQ(310, (*after)(6,9,0));
+   }
+   // Mask in values
+   TEST_F(TestCalibratorMask, 10x10_mask_in) {
+      FileArome from("testing/files/10x10.nc");
+      ParameterFileText par(Options("file=testing/files/mask0.txt spatial=1"));
+      CalibratorMask cal = CalibratorMask(Variable::T, Options("mask=1"));
+
+      cal.calibrate(from, &par);
+      FieldPtr after = from.getField(Variable::T, 0);
+      ASSERT_EQ(10, after->getNumLat());
+      ASSERT_EQ(10, after->getNumLon());
+      ASSERT_EQ(1,  after->getNumEns());
+
+      EXPECT_FLOAT_EQ(Util::MV, (*after)(5,2,0));
+      EXPECT_FLOAT_EQ(302, (*after)(3,5,0));
+      EXPECT_FLOAT_EQ(302, (*after)(3,3,0));
+      EXPECT_FLOAT_EQ(302, (*after)(2,5,0));
+      EXPECT_FLOAT_EQ(302, (*after)(4,9,0));
+      EXPECT_FLOAT_EQ(Util::MV, (*after)(2,3,0));
+      EXPECT_FLOAT_EQ(Util::MV, (*after)(6,9,0));
+   }
+   // Missing parameter file
+   TEST_F(TestCalibratorMask, invalid2) {
+      ::testing::FLAGS_gtest_death_test_style = "threadsafe";
+	  FileArome from("testing/files/10x10.nc");
+      Util::setShowError(false);
+      CalibratorMask calibrator(Variable::T, Options());
+      EXPECT_DEATH(calibrator.calibrate(from, NULL), ".*");
+   }
+   TEST_F(TestCalibratorMask, description) {
+      CalibratorMask::description();
+   }
+}
+int main(int argc, char **argv) {
+     ::testing::InitGoogleTest(&argc, argv);
+       return RUN_ALL_TESTS();
+}

src/Version.h

@@ -1,3 +1,3 @@
 #ifndef GRIDPP_VERSION
-#define GRIDPP_VERSION "0.2.6"
+#define GRIDPP_VERSION "0.2.7"
 #endif

testing/files/mask0.txt

@@ -0,0 +1,2 @@
+0 3 5 0 223000
+0 4 6 0 336000