Ticket 779: Correction of some issues, see trac for more info jenkins-beamb-66
authorulf.nordh <a002011@c21601.ad.smhi.se>
Wed, 5 Jun 2019 11:16:38 +0000 (13:16 +0200)
committerulf.nordh <a002011@c21601.ad.smhi.se>
Wed, 5 Jun 2019 11:16:38 +0000 (13:16 +0200)
bin/beamb
lib/beamblockage.c
pybeamb/beamb_quality_plugin.py

index 541fab5..98703a7 100755 (executable)
--- a/bin/beamb
+++ b/bin/beamb
@@ -42,7 +42,7 @@ def process_scan(scan, options):
         restored = _beamblockage.restore(scan, result, options.quantity, options.restore)
         scan.addQualityField(result)
     else:
-        print "Ignoring scan with elevation angle %3.1f degrees" % (scan.elangle*rd)
+        print("Ignoring scan with elevation angle %3.1f degrees"%(scan.elangle*rd))
     return scan
 
 
@@ -87,7 +87,7 @@ if __name__ == "__main__":
                       help="Specifies the limit of the Gaussian approximation of the beamwidth (in dB) to apply when analyzing the degree of beam blockage. NOTE that this does not have to be the radar's half-power beamwidth. Defaults to -6.0 dB.")
 
     parser.add_option("-r", "--restore", dest="restore", default=0.7, type="float",
-                      help="Specifies the upper threshold on how much beam blockage is accepted to be corrected. Sectors blocked more than this value are blocked out by assigning them the 'nodata' value. Defaults to 0.7 (70 %). Values above 1 will correct all data and mask none.")
+                      help="Specifies the upper threshold on how much beam blockage is accepted to be corrected. Sectors blocked more than this value are blocked out by assigning them the 'nodata' value. Defaults to 0.7 (70 %). Values needs to be between 0.0 and 1.0, otherwise the code will fail to restore. To avoid having masking with nodata, the threshold shall be set to 1.0")
     
     parser.add_option("-q", "--quantity", dest="quantity", default="DBZH",
                       help="Specifies the quantity to work with, default 'DBZH'.")
@@ -101,7 +101,7 @@ if __name__ == "__main__":
         rio = _raveio.open(options.infile)
 
         if rio.objectType not in (_raveio.Rave_ObjectType_PVOL, _raveio.Rave_ObjectType_SCAN):
-            print "Input file must be either polar scan or volume. Exiting ..."
+            print("Input file must be either polar scan or volume. Exiting ...")
             sys.exit(1)
 
         elif rio.objectType == _raveio.Rave_ObjectType_PVOL:
index f4bfd54..7e8fc81 100644 (file)
@@ -23,6 +23,9 @@ along with beamb.  If not, see <http://www.gnu.org/licenses/>.
  *
  * @author Anders Henja (SMHI, refactored to work together with rave)
  * @date 2011-11-10
+ *
+ * @author Ulf E. Nordh (SMHI, removal of bugs in fkns computeGroundRange and restore)
+ * @date 2019-05-22 
  */
 #include "beamblockage.h"
 #include "beamblockagemap.h"
@@ -126,6 +129,7 @@ static double* BeamBlockageInternal_computeGroundRange(BeamBlockage_t* self, Pol
   navigator = PolarScan_getNavigator(scan);
   nbins = PolarScan_getNbins(scan);
   rscale = PolarScan_getRscale(scan);
+  elangle = PolarScan_getElangle(scan);
 
   ranges = RAVE_MALLOC(sizeof(double) * nbins);
   if (ranges == NULL) {
@@ -572,11 +576,13 @@ RaveField_t* BeamBlockage_getBlockage(BeamBlockage_t* self, PolarScan_t* scan, d
       }
       bbval = -1.0/2.0 * sqrt(M_PI * c) * (erf((elangle - elBlock)/sqrt(c)) - erf(elLim/sqrt(c)))/bb_tot;
 
+      /* Discard non-physical values for blockage */
       if (bbval < 0.0) {
         bbval = 0.0;
       } else if (bbval > 1.0) {
         bbval = 1.0;
       }
+
       /* ODIM's rule for representing quality is that 0=lowest, 1=highest quality. Therefore invert. */
       bbval = ((1.0-bbval) - offset) / gain;
       RaveField_setValue(field, bi, ri, bbval);
@@ -607,8 +613,8 @@ int BeamBlockage_restore(PolarScan_t* scan, RaveField_t* blockage, const char* q
   PolarScanParam_t* parameter = NULL;
   double gain, offset, nodata;
   long ri, bi, nrays, nbins;
-  double bbgain, bboffset;
-  double iv, ov, bbraw, bbpercent, bbdbz, dbz;
+  double bbgain, bboffset, bbraw;
+  double dbz_uncorr, dbz_corr, dbz_corr_raw, bbpercent, bb_corr_db, bb_corr_lin;
   RaveValueType rvt;
 
   if (scan == NULL || blockage == NULL) {
@@ -642,26 +648,50 @@ int BeamBlockage_restore(PolarScan_t* scan, RaveField_t* blockage, const char* q
     goto done;
   }
 
+  /* Exits if the user failed to set the threshold between 0.0 and 1.0 */
+  if (threshold < 0.0 || threshold > 1.0) {
+    RAVE_ERROR0("A blockage threshold smaller than 0.0 or larger than 1.0 is set in the calling script, correct and retry");
+    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) || (rvt == RaveValueType_UNDETECT)) {
+
+      /* The function below does two things, 1) checks the data type of the matrix element, and 2) converts from raw to dbz */
+      /* Returns in variable rvt and the variable dbz_uncorr */
+      rvt = PolarScanParam_getConvertedValue(parameter, bi, ri, &dbz_uncorr);
+
+      if ((rvt == RaveValueType_DATA) || (rvt == RaveValueType_UNDETECT)) { /* Ignore nodata matrix elements */
         RaveField_getValue(blockage, bi, ri, &bbraw);
+
         /* ODIM's rule for representing quality is that 0=lowest, 1=highest quality. Therefore revert. */
         bbpercent = 1.0 - (bbgain * bbraw + bboffset);
+
+        /* Discard non-physical values */
         if (bbpercent < 0.0 || bbpercent > 1.0) {
           RAVE_ERROR0("beamb values are out of bounds, check scaling");
           goto done;
         }
-        if ((bbpercent < threshold) && (rvt == RaveValueType_DATA)) {
-          bbdbz = 10 * (log10(1.0/(pow(1.0-bbpercent,2)))); /* Two-way correction */
-          dbz = iv + bbdbz;  /* This is the corrected reflectivity */
-          ov = round((dbz - offset) / gain);
-
-/*           if (ov > nodata - 1) */ /* Risky to use and not to use... */
-/*             ov = nodata - 1; */
-          PolarScanParam_setValue(parameter, bi, ri, ov);
-        } else if (bbpercent >= threshold) {
+
+        /* Adjust the reflectivity IF we have blockage and IF it is smaller than the selected threshold AND we are dealing with data */
+        if ((bbpercent > 0.0) && (bbpercent <= threshold) && (bbpercent != 1.0) && (rvt == RaveValueType_DATA)) {
+
+          bb_corr_lin = 1.0 / (pow(1.0-bbpercent,2)); /* Two-way blockage multiplicative correction in linear representation */
+          bb_corr_db = 10.0 * log10(bb_corr_lin); /* Two-way blockage multiplicative correction in dB */
+
+          dbz_corr = dbz_uncorr + bb_corr_db;  /* In the "dB-regime" we use addition for the multiplicative correction */
+          dbz_corr_raw = round((dbz_corr - offset) / gain); /* Converting to raw format */
+
+          /* if (dbz_corr_raw > nodata - 1)  Brute force used to in worst case keep the corrected data within 8-bit,perhaps risky to use 
+          {
+             dbz_corr_raw = nodata - 1;
+          } */
+
+          PolarScanParam_setValue(parameter, bi, ri, dbz_corr_raw);
+
+        /* If the blockage is larger than the selected threshold, mask with nodata, this is valid for both data and undetect */
+        /* By doing like this we give adjacent radars the opportunity to fill in these pixels */
+        } else if (bbpercent > threshold) {
           PolarScanParam_setValue(parameter, bi, ri, nodata);  /* Uncorrectable */
         }
       }
index c6e08e3..0d310fe 100644 (file)
@@ -40,7 +40,14 @@ logger = rave_pgf_logger.create_logger()
 # The limit of the Gaussian approximation of main lobe
 #
 BEAMBLOCKAGE_DBLIMIT=-6.0
-BEAMBLOCKAGE_BBLIMIT= 1.1
+
+##
+# The threshold for correcting beamblockage, must be between 0.0 and 1.0, if set to 1.0,
+# no masking with nodata will be done i.e correction will be done for all blockage up to 1.0
+# with 1.0 excluded due to that the underlying C-code cannot deal with 1.0 (division by zero). An attempt to set
+# the parameter smaller than 0.0 or larger than 1.0 will cause the restore function to abort.
+#
+BEAMBLOCKAGE_BBLIMIT= 1.0 
 
 ##
 # The beam blockage quality plugin
@@ -64,8 +71,7 @@ class beamb_quality_plugin(rave_quality_plugin):
   _dblimit = BEAMBLOCKAGE_DBLIMIT
 
   ##
-  # The default percent beam blockage (divided by 100). Defaults to BEAMBLOCKAGE_BBLIMIT which 
-  # is set to 110% so that no radar data will be masked to NODATA.
+  # The default percent beam blockage (divided by 100). Defaults to BEAMBLOCKAGE_BBLIMIT.
   #
   _bblimit = BEAMBLOCKAGE_BBLIMIT