Ticket 741: Add restore functionality to the beam blockage module hudson-beamb-27-SUCCESS
authorAnders Henja <anders@baltrad.eu>
Wed, 4 Jan 2012 15:27:37 +0000 (16:27 +0100)
committerAnders Henja <anders@baltrad.eu>
Wed, 4 Jan 2012 15:27:37 +0000 (16:27 +0100)
lib/beamblockage.c
lib/beamblockage.h
pybeamb/pybeamblockage.c
test/pytest/PyBeamBlockageTest.py

index eb4cf56..caf58d4 100644 (file)
@@ -44,6 +44,12 @@ struct _BeamBlockage_t {
   int rewritecache;         /**< if cache should be recreated */
 };
 
+/**
+ * Converts a radian to a degree
+ * @param[in] rad - input value expressed in radians
+ */
+#define RAD2DEG(rad) (rad*180.0/M_PI)
+
 /*@{ Private functions */
 /**
  * Constructor.
@@ -253,6 +259,50 @@ done:
   return result;
 }
 
+static int BeamBlockageInternal_getMetaInformation(RaveField_t* field, double* gain, double* offset)
+{
+  int result = 0;
+  RaveAttribute_t* attribute = NULL;
+  char* svalue = NULL;
+
+  RAVE_ASSERT((field != NULL), "field == NULL");
+  RAVE_ASSERT((gain != NULL), "gain == NULL");
+  RAVE_ASSERT((offset != NULL), "offset == NULL");
+
+  attribute = RaveField_getAttribute(field, "how/task");
+  if (attribute == NULL || !RaveAttribute_getString(attribute, &svalue)) {
+    RAVE_ERROR0("how/task missing");
+    goto done;
+  }
+  if (svalue == NULL || strcmp("se.smhi.detector.beamblockage", svalue) != 0) {
+    if (svalue == NULL) {
+      RAVE_ERROR0("how/task == NULL");
+    } else {
+      RAVE_ERROR1("how/task = %s", svalue);
+    }
+    goto done;
+  }
+  RAVE_OBJECT_RELEASE(attribute);
+
+  attribute = RaveField_getAttribute(field, "what/gain");
+  if (attribute == NULL || !RaveAttribute_getDouble(attribute, gain)) {
+    RAVE_ERROR0("Missing what/gain");
+    goto done;
+  }
+  RAVE_OBJECT_RELEASE(attribute);
+
+  attribute = RaveField_getAttribute(field, "what/offset");
+  if (attribute == NULL || !RaveAttribute_getDouble(attribute, offset)) {
+    RAVE_ERROR0("Missing what/offset");
+    goto done;
+  }
+
+  result = 1;
+done:
+  RAVE_OBJECT_RELEASE(attribute);
+  return result;
+}
+
 /**
  * Returns a cached file matching the given scan if there is one.
  * @param[in] self - self
@@ -355,16 +405,6 @@ done:
   return result;
 }
 
-/**
- * Convert radians to degrees
- * @param[in] rad - input value expressed in radians
- * @returns value converted to degrees
- */
-double rad2deg(double rad)
-{
-    return rad*180./M_PI;
-}
-
 /*@} End of Private functions */
 
 /*@{ Interface functions */
@@ -501,7 +541,7 @@ RaveField_t* BeamBlockage_getBlockage(BeamBlockage_t* self, PolarScan_t* scan, d
     for (bi = 0; bi < nbins; bi++) {
       double v = 0.0;
       BBTopography_getValue(topo, bi, ri, &v);
-      phi[bi] = rad2deg(asin((((v+R)*(v+R)) - (groundRange[bi]*groundRange[bi]) - ((R+height)*(R+height))) / (2*groundRange[bi]*(R+height))));
+      phi[bi] = RAD2DEG(asin((((v+R)*(v+R)) - (groundRange[bi]*groundRange[bi]) - ((R+height)*(R+height))) / (2*groundRange[bi]*(R+height))));
     }
     BeamBlockageInternal_cummax(phi, nbins);
 
@@ -544,6 +584,78 @@ done:
   return result;
 }
 
+int BeamBlockage_restore(PolarScan_t* scan, RaveField_t* blockage, const char* quantity, double threshold)
+{
+  int result = 0;
+  PolarScanParam_t* parameter = NULL;
+  double gain, offset, nodata;
+  long ri, bi, nrays, nbins;
+  double bbgain, bboffset;
+  double iv, ov, bbraw, bbpercent, bbdbz, dbz;
+  RaveValueType rvt;
+
+  if (scan == NULL || blockage == NULL) {
+    RAVE_ERROR0("Need to provide both scan and field containing blockage.");
+    goto done;
+  }
+
+  if (!BeamBlockageInternal_getMetaInformation(blockage, &bbgain, &bboffset)) {
+    RAVE_ERROR0("Could not get meta information from blockage field.");
+    goto done;
+  }
+
+  if (quantity == NULL) {
+    parameter = PolarScan_getParameter(scan, "DBZH");
+  } else {
+    parameter = PolarScan_getParameter(scan, quantity);
+  }
+  if (parameter == NULL) {
+    RAVE_ERROR1("No parameter with quantity %s in scan", quantity);
+    goto done;
+  }
+
+  gain = PolarScanParam_getGain(parameter);
+  offset = PolarScanParam_getOffset(parameter);
+  nodata = PolarScanParam_getNodata(parameter);
+  nrays = PolarScanParam_getNrays(parameter);
+  nbins = PolarScanParam_getNbins(parameter);
+
+  if (nrays != RaveField_getYsize(blockage) || nbins != RaveField_getXsize(blockage)) {
+    RAVE_ERROR0("field and scan dimensions must be the same");
+    goto done;
+  }
+
+  for (ri = 0; ri < nrays; ri++) {
+    for (bi = 0; bi < nbins; bi++) {
+      rvt = PolarScanParam_getConvertedValue(parameter, bi, ri, &iv);
+      if (rvt == RaveValueType_DATA) {
+        RaveField_getValue(blockage, bi, ri, &bbraw);
+        bbpercent = bbgain * bbraw + bboffset;
+        if (bbpercent < 0.0 || bbpercent > 1.0) {
+          RAVE_ERROR0("beamb values are out of bounds, check scaling");
+          goto done;
+        }
+        if (bbpercent < threshold) {
+          bbdbz = (10*log(-1.0/(bbpercent-1)))/log(10.);
+          dbz = iv + bbdbz;  /* This is the corrected reflectivity */
+          ov = (dbz - offset) / gain;
+
+          if (ov > nodata - 1)
+            ov = nodata - 1;
+          PolarScanParam_setValue(parameter, bi, ri, ov);
+        } else {
+          PolarScanParam_setValue(parameter, bi, ri, nodata);  /* Uncorrectable */
+        }
+      }
+    }
+  }
+
+  result = 1;
+done:
+  RAVE_OBJECT_RELEASE(parameter);
+  return result;
+}
+
 /*@} End of Interface functions */
 
 RaveCoreObjectType BeamBlockage_TYPE = {
index f923129..21d95f4 100644 (file)
@@ -97,4 +97,14 @@ int BeamBlockage_getRewriteCache(BeamBlockage_t* self);
  */
 RaveField_t* BeamBlockage_getBlockage(BeamBlockage_t* self, PolarScan_t* scan, double dBlim);
 
+/**
+ * When you have retrieved the beam blockage field you can restore the specified parameter
+ * for the scan.
+ * @param[in] scan - the scan that was provided to the getBlockage function
+ * @param[in] blockage - the result from the call to getBlockage
+ * @param[in] quantity - the parameter to be restored. If NULL, it defaults to DBZH
+ * @param[in] threshold - the percentage threshold
+ */
+int BeamBlockage_restore(PolarScan_t* scan, RaveField_t* blockage, const char* quantity, double threshold);
+
 #endif /* BEAMBLOCKAGE_H */
index a9edd09..e9cf653 100644 (file)
@@ -146,6 +146,35 @@ static PyObject* _pybeamblockage_new(PyObject* self, PyObject* args)
 }
 
 /**
+ * Restores the provided scan with the beam blockage field.
+ * @param[in] self - this instance
+ * @param[in] args - (OOsf), (scan, field, quantity, threshold value)
+ * @returns None
+ */
+static PyObject* _pybeamblockage_restore(PyObject* self, PyObject* args)
+{
+  PyObject *o1 = NULL, *o2 = NULL;
+  char* quantity = NULL;
+  double threshold = 0.0;
+
+  if (!PyArg_ParseTuple(args, "OOsf", &o1, &o2, &quantity, &threshold)) {
+    return NULL;
+  }
+
+  if (!PyPolarScan_Check(o1)) {
+    raiseException_returnNULL(PyExc_TypeError, "First argument should be a PolarScan");
+  }
+  if (!PyRaveField_Check(o2)) {
+    raiseException_returnNULL(PyExc_TypeError, "Second argument should be a PolarScan");
+  }
+
+  if (!BeamBlockage_restore(((PyPolarScan*)o1)->scan, ((PyRaveField*)o2)->field, quantity, threshold)) {
+    raiseException_returnNULL(PyExc_RuntimeError, "Failed to restore scan");
+  }
+  Py_RETURN_NONE;
+}
+
+/**
  * Returns the blockage for the provided scan given gaussian limit.
  * @param[in] self - self
  * @param[in] args - the arguments (PyPolarScan, double (Limit of Gaussian approximation of main lobe))
@@ -299,6 +328,7 @@ PyTypeObject PyBeamBlockage_Type =
 /*@{ Module setup */
 static PyMethodDef functions[] = {
   {"new", (PyCFunction)_pybeamblockage_new, 1},
+  {"restore", (PyCFunction)_pybeamblockage_restore, 1},
   {NULL,NULL} /*Sentinel*/
 };
 
index 8f71451..1f12e29 100644 (file)
@@ -29,6 +29,7 @@ import _raveio
 import _beamblockage
 import os, string
 import _rave
+import _ravefield
 import numpy
 
 class PyBeamBlockageTest(unittest.TestCase):
@@ -165,6 +166,69 @@ class PyBeamBlockageTest(unittest.TestCase):
     self.assertTrue(os.path.isfile(self.CACHEFILE_3))
     c2stat = os.stat(self.CACHEFILE_2)
 
+  def test_restore(self):
+    a = _beamblockage.new()
+    a.topo30dir="../../data/gtopo30"
+    a.cachedir="/tmp"
+    a.rewritecache=True
+    scan = _raveio.open(self.FIXTURE_2).object
+    blockage = a.getBlockage(scan, -20.0);
+
+    restored = _beamblockage.restore(scan, blockage, "DBZH", 0.5)
+
+  def test_restore_proper_field(self):
+    scan = _raveio.open(self.FIXTURE_2).object
+
+    field = _ravefield.new()
+    field.setData(numpy.zeros((scan.nrays, scan.nbins), numpy.uint8))
+    field.addAttribute("how/task", "se.smhi.detector.beamblockage")
+    field.addAttribute("what/gain", 1.0)
+    field.addAttribute("what/offset", 0.0)
+        
+    _beamblockage.restore(scan, field, "DBZH", 0.5)
+    
+  def test_restore_missing_howtask(self):
+    scan = _raveio.open(self.FIXTURE_2).object
+
+    field = _ravefield.new()
+    field.setData(numpy.zeros((scan.nrays, scan.nbins), numpy.uint8))
+    field.addAttribute("what/gain", 1.0)
+    field.addAttribute("what/offset", 0.0)
+
+    try:        
+      _beamblockage.restore(scan, field, "DBZH", 0.5)
+      self.fail("Expected RuntimeError due to missing how/task")
+    except RuntimeError:
+      pass
+    
+  def test_restore_missing_whatgain(self):
+    scan = _raveio.open(self.FIXTURE_2).object
+
+    field = _ravefield.new()
+    field.setData(numpy.zeros((scan.nrays, scan.nbins), numpy.uint8))
+    field.addAttribute("how/task", "se.smhi.detector.beamblockage")
+    field.addAttribute("what/offset", 0.0)
+
+    try:        
+      _beamblockage.restore(scan, field, "DBZH", 0.5)
+      self.fail("Expected RuntimeError due to missing how/gain")
+    except RuntimeError:
+      pass
+
+  def test_restore_missing_whatoffset(self):
+    scan = _raveio.open(self.FIXTURE_2).object
+
+    field = _ravefield.new()
+    field.setData(numpy.zeros((scan.nrays, scan.nbins), numpy.uint8))
+    field.addAttribute("how/task", "se.smhi.detector.beamblockage")
+    field.addAttribute("what/gain", 1.0)
+
+    try:        
+      _beamblockage.restore(scan, field, "DBZH", 0.5)
+      self.fail("Expected RuntimeError due to missing how/offset")
+    except RuntimeError:
+      pass
+    
 if __name__ == "__main__":
   #import sys;sys.argv = ['', 'Test.testName']
   unittest.main()
\ No newline at end of file