lib/wrwp.c

   1 /* --------------------------------------------------------------------
   2 Copyright (C) 2011 Swedish Meteorological and Hydrological Institute, SMHI
   3
   4 This is free software: you can redistribute it and/or modify
   5 it under the terms of the GNU Lesser General Public License as published by
   6 the Free Software Foundation, either version 3 of the License, or
   7 (at your option) any later version.
   8
   9 This software is distributed in the hope that it will be useful,
  10 but WITHOUT ANY WARRANTY; without even the implied warranty of
  11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12 GNU Lesser General Public License for more details.
  13
  14 You should have received a copy of the GNU Lesser General Public License
  15 along with HLHDF.  If not, see <http://www.gnu.org/licenses/>.
  16 ------------------------------------------------------------------------*/
  17
  18 /** Function for deriving weather radar wind and reflectivity profiles
  19  *
  20  * @file
  21  * @author Gunther Haase, SMHI
  22  * @date 2013-02-06
  23  */
  24
  25 #include "wrwp.h"
  26 #include "rave_debug.h"
  27 #include "rave_alloc.h"
  28 #include <string.h>
  29
  30 /**
  31  * Represents one wrwp generator
  32  */
  33 struct _Wrwp_t {
  34   RAVE_OBJECT_HEAD /** Always on top */
  35   int dz; /**< Height interval for deriving a profile [m] */
  36   int hmax; /**< Maximum height of the profile [m] */
  37   int dmin; /**< Minimum distance for deriving a profile [m] */
  38   int dmax; /**< Maximum distance for deriving a profile [m]*/
  39   double emin; /**< Minimum elevation angle [deg] */
  40   double vmin; /**< Radial velocity threshold [m/s] */
  41 };
  42
  43 /*@{ Private functions */
  44 /**
  45  * Constructor
  46  */
  47 static int Wrwp_constructor(RaveCoreObject* obj)
  48 {
  49   Wrwp_t* wrwp = (Wrwp_t*)obj;
  50   wrwp->hmax = HMAX;
  51   wrwp->dz = DZ;
  52   wrwp->dmin = DMIN;
  53   wrwp->dmax = DMAX;
  54   wrwp->emin = EMIN;
  55   wrwp->vmin = VMIN;
  56   return 1;
  57 }
  58
  59 /**
  60  * Destructor
  61  */
  62 static void Wrwp_destructor(RaveCoreObject* obj)
  63 {
  64 }
  65
  66 static int WrwpInternal_findAndAddAttribute(VerticalProfile_t* vp, PolarVolume_t* pvol, const char* name)
  67 {
  68   int nscans = PolarVolume_getNumberOfScans(pvol);
  69   int i = 0;
  70   int found = 0;
  71   for (i = 0; i < nscans && found == 0; i++) {
  72     PolarScan_t* scan = PolarVolume_getScan(pvol, i);
  73     if (scan != NULL && PolarScan_hasAttribute(scan, name)) {
  74       RaveAttribute_t* attr = PolarScan_getAttribute(scan, name);
  75       VerticalProfile_addAttribute(vp, attr);
  76       found = 1;
  77       RAVE_OBJECT_RELEASE(attr);
  78     }
  79     RAVE_OBJECT_RELEASE(scan);
  80   }
  81   return found;
  82 }
  83
  84 static int WrwpInternal_addIntAttribute(VerticalProfile_t* vp, const char* name, int value)
  85 {
  86   RaveAttribute_t* attr = RaveAttributeHelp_createLong(name, value);
  87   int result = 0;
  88   if (attr != NULL) {
  89     result = VerticalProfile_addAttribute(vp, attr);
  90   }
  91   RAVE_OBJECT_RELEASE(attr);
  92   return result;
  93 }
  94
  95 /*@} End of Private functions */
  96
  97 /*@{ Interface functions */
  98 void Wrwp_setDZ(Wrwp_t* self, int dz)
  99 {
 100   RAVE_ASSERT((self != NULL), "self == NULL");
 101   self->dz = dz;
 102 }
 103
 104 int Wrwp_getDZ(Wrwp_t* self)
 105 {
 106   RAVE_ASSERT((self != NULL), "self == NULL");
 107   return self->dz;
 108 }
 109
 110 void Wrwp_setHMAX(Wrwp_t* self, int hmax)
 111 {
 112   RAVE_ASSERT((self != NULL), "self == NULL");
 113   self->hmax = hmax;
 114 }
 115
 116 int Wrwp_getHMAX(Wrwp_t* self)
 117 {
 118   RAVE_ASSERT((self != NULL), "self == NULL");
 119   return self->hmax;
 120 }
 121
 122 void Wrwp_setDMIN(Wrwp_t* self, int dmin)
 123 {
 124   RAVE_ASSERT((self != NULL), "self == NULL");
 125   self->dmin = dmin;
 126 }
 127
 128 int Wrwp_getDMIN(Wrwp_t* self)
 129 {
 130   RAVE_ASSERT((self != NULL), "self == NULL");
 131   return self->dmin;
 132 }
 133
 134 void Wrwp_setDMAX(Wrwp_t* self, int dmax)
 135 {
 136   RAVE_ASSERT((self != NULL), "self == NULL");
 137   self->dmax = dmax;
 138 }
 139
 140 int Wrwp_getDMAX(Wrwp_t* self)
 141 {
 142   RAVE_ASSERT((self != NULL), "self == NULL");
 143   return self->dmax;
 144 }
 145
 146 void Wrwp_setEMIN(Wrwp_t* self, double emin)
 147 {
 148   RAVE_ASSERT((self != NULL), "self == NULL");
 149   self->emin = emin;
 150 }
 151
 152 double Wrwp_getEMIN(Wrwp_t* self)
 153 {
 154   RAVE_ASSERT((self != NULL), "self == NULL");
 155   return self->emin;
 156 }
 157
 158 void Wrwp_setVMIN(Wrwp_t* self, double vmin)
 159 {
 160   RAVE_ASSERT((self != NULL), "self == NULL");
 161   self->vmin = vmin;
 162 }
 163
 164 double Wrwp_getVMIN(Wrwp_t* self)
 165 {
 166   RAVE_ASSERT((self != NULL), "self == NULL");
 167   return self->vmin;
 168 }
 169
 170 VerticalProfile_t* Wrwp_generate(Wrwp_t* self, PolarVolume_t* inobj) {
 171   VerticalProfile_t* result = NULL;
 172         PolarScan_t* scan = NULL;
 173         PolarScanParam_t* vrad = NULL;
 174         PolarScanParam_t* dbz = NULL;
 175         PolarNavigator_t* polnav = NULL;
 176         int nrhs = NRHS, lda = LDA, ldb = LDB;
 177         int nscans = 0, nbins = 0, nrays, nv, nz, i, iz, is, ib, ir;
 178         double rscale, elangle, gain, offset, nodata, undetect, val;
 179         double d, h;
 180         double alpha, beta, gamma, vvel, vdir, vstd, zsum, zmean, zstd;
 181         int ysize = 0, yindex = 0;
 182
 183         RaveField_t *ff_field = NULL, *ff_dev_field = NULL, *dd_field = NULL;
 184         RaveField_t *dbzh_field = NULL, *dbzh_dev_field = NULL, *nz_field = NULL, *nv_field = NULL;
 185
 186         RAVE_ASSERT((self != NULL), "self == NULL");
 187
 188         ff_field = RAVE_OBJECT_NEW(&RaveField_TYPE);
 189   ff_dev_field = RAVE_OBJECT_NEW(&RaveField_TYPE);
 190   dd_field = RAVE_OBJECT_NEW(&RaveField_TYPE);
 191   dbzh_field = RAVE_OBJECT_NEW(&RaveField_TYPE);
 192   dbzh_dev_field = RAVE_OBJECT_NEW(&RaveField_TYPE);
 193   nz_field = RAVE_OBJECT_NEW(&RaveField_TYPE);
 194   nv_field = RAVE_OBJECT_NEW(&RaveField_TYPE);
 195
 196   if (ff_field == NULL || ff_dev_field == NULL || dd_field == NULL ||
 197       dbzh_field == NULL || dbzh_dev_field == NULL || nz_field == NULL || nv_field == NULL) {
 198     RAVE_ERROR0("Failed to allocate memory for the resulting vp fields");
 199     goto done;
 200   }
 201   ysize = self->hmax / self->dz;
 202   if (!RaveField_createData(ff_field, 1, ysize, RaveDataType_DOUBLE) ||
 203       !RaveField_createData(ff_dev_field, 1, ysize, RaveDataType_DOUBLE) ||
 204       !RaveField_createData(dd_field, 1, ysize, RaveDataType_DOUBLE) ||
 205       !RaveField_createData(dbzh_field, 1, ysize, RaveDataType_DOUBLE) ||
 206       !RaveField_createData(dbzh_dev_field, 1, ysize, RaveDataType_DOUBLE) ||
 207       !RaveField_createData(nz_field, 1, ysize, RaveDataType_INT) ||
 208       !RaveField_createData(nv_field, 1, ysize, RaveDataType_INT)) {
 209     RAVE_ERROR0("Failed to allocate arrays for the resulting vp fields");
 210     goto done;
 211   }
 212
 213         polnav = RAVE_OBJECT_NEW (&PolarNavigator_TYPE);
 214         PolarNavigator_setLat0(polnav, PolarVolume_getLatitude(inobj));
 215         PolarNavigator_setLon0(polnav, PolarVolume_getLongitude(inobj));
 216         PolarNavigator_setAlt0(polnav, PolarVolume_getHeight(inobj));
 217
 218         nscans = PolarVolume_getNumberOfScans (inobj);
 219
 220         // We use yindex for filling in the arrays even though we loop to hmax...
 221         yindex = 0;
 222
 223         // loop over atmospheric layers
 224   for (iz = 0; iz < self->hmax; iz += self->dz) {
 225     /* allocate memory and initialize with zeros */
 226     double *A = RAVE_CALLOC((size_t)(NOR*NOC), sizeof (double));
 227     double *b = RAVE_CALLOC((size_t)(NOR), sizeof (double));
 228     double *v = RAVE_CALLOC((size_t)(NOR), sizeof (double));
 229     double *z = RAVE_CALLOC((size_t)(NOR), sizeof (double));
 230     double *az = RAVE_CALLOC((size_t)(NOR), sizeof (double));
 231
 232     vdir = -9999.0;
 233     vvel = -9999.0;
 234     vstd = 0.0;
 235     zsum = 0.0;
 236     zmean = -9999.0;
 237     zstd = 0.0;
 238     nv = 0;
 239     nz = 0;
 240
 241     // loop over scans
 242     for (is = 0; is < nscans; is++) {
 243       scan = PolarVolume_getScan(inobj, is);
 244       nbins = PolarScan_getNbins(scan);
 245       nrays = PolarScan_getNrays(scan);
 246       rscale = PolarScan_getRscale(scan);
 247       elangle = PolarScan_getElangle(scan);
 248
 249       // radial wind scans
 250       if (PolarScan_hasParameter(scan, "VRAD")) {
 251         vrad = PolarScan_getParameter(scan, "VRAD");
 252         gain = PolarScanParam_getGain(vrad);
 253         offset = PolarScanParam_getOffset(vrad);
 254         nodata = PolarScanParam_getNodata(vrad);
 255         undetect = PolarScanParam_getUndetect(vrad);
 256
 257         for (ir = 0; ir < nrays; ir++) {
 258           for (ib = 0; ib < nbins; ib++) {
 259             PolarNavigator_reToDh(polnav, (ib+0.5)*rscale, elangle, &d, &h);
 260             PolarScanParam_getValue(vrad, ib, ir, &val);
 261
 262             if ((h >= iz) &&
 263                 (h < iz + self->dz) &&
 264                 (d >= self->dmin) &&
 265                 (d <= self->dmax) &&
 266                 (elangle * RAD2DEG >= self->emin) &&
 267                 (val != nodata) &&
 268                 (val != undetect) &&
 269                 (abs(offset + gain * val) >= self->vmin)) {
 270               *(v+nv) = offset+gain*val;
 271               *(az+nv) = 360./nrays*ir*DEG2RAD;
 272               *(A+nv*NOC) = sin(*(az+nv));
 273               *(A+nv*NOC+1) = cos(*(az+nv));
 274               *(A+nv*NOC+2) = 1;
 275               *(b+nv) = *(v+nv);
 276               nv = nv+1;
 277             }
 278           }
 279         }
 280         RAVE_OBJECT_RELEASE(vrad);
 281       }
 282
 283       // reflectivity scans
 284       if (PolarScan_hasParameter(scan, "DBZH")) {
 285         dbz = PolarScan_getParameter(scan, "DBZH");
 286         gain = PolarScanParam_getGain(dbz);
 287         offset = PolarScanParam_getOffset(dbz);
 288         nodata = PolarScanParam_getNodata(dbz);
 289         undetect = PolarScanParam_getUndetect(dbz);
 290
 291         for (ir = 0; ir < nrays; ir++) {
 292           for (ib = 0; ib < nbins; ib++) {
 293             PolarNavigator_reToDh (polnav, (ib+0.5)*rscale, elangle, &d, &h);
 294             PolarScanParam_getValue (dbz, ib, ir, &val);
 295             if ((h >= iz) &&
 296                 (h < iz+self->dz) &&
 297                 (d >= self->dmin) &&
 298                 (d <= self->dmax) &&
 299                 (elangle*RAD2DEG >= self->emin) &&
 300                 (val != nodata) &&
 301                 (val != undetect)) {
 302               *(z+nz) = dBZ2Z(offset+gain*val);
 303               zsum = zsum + *(z+nz);
 304               nz = nz+1;
 305             }
 306           }
 307         }
 308         RAVE_OBJECT_RELEASE(dbz);
 309       }
 310       RAVE_OBJECT_RELEASE(scan);
 311     }
 312
 313     // radial wind calculations
 314     if (nv>3) {
 315       //***************************************************************
 316       // fitting: y = gamma+alpha*sin(x+beta)                         *
 317       // alpha -> amplitude                                           *
 318       // beta -> phase shift                                          *
 319       // gamma -> consider an y-shift due to the terminal velocity of *
 320       //          falling rain drops                                  *
 321       //***************************************************************
 322
 323       /* QR decomposition */
 324       /*info = */LAPACKE_dgels(LAPACK_ROW_MAJOR, 'N', NOR, NOC, nrhs, A, lda, b, ldb);
 325
 326       /* parameter of the wind model */
 327       alpha = sqrt(pow(*(b),2) + pow(*(b+1),2));
 328       beta = atan2(*(b+1), *b);
 329       gamma = *(b+2);
 330
 331       /* wind velocity */
 332       vvel = alpha;
 333
 334       /* wind direction */
 335       vdir = 0;
 336       if (alpha < 0) {
 337         vdir = (PI/2-beta) * RAD2DEG;
 338       } else if (alpha > 0) {
 339         vdir = (3*PI/2-beta) * RAD2DEG;
 340       }
 341       if (vdir < 0) {
 342         vdir = vdir + 360;
 343       } else if (vdir > 360) {
 344         vdir = vdir - 360;
 345       }
 346
 347       /* standard deviation of the wind velocity */
 348       for (i=0; i<nv; i++) {
 349         vstd = vstd + pow (*(v+i) - (gamma+alpha*sin (*(az+i)+beta)),2);
 350       }
 351       vstd = sqrt (vstd/(nv-1));
 352     }
 353
 354     // reflectivity calculations
 355     if (nz > 1) {
 356       /* standard deviation of reflectivity */
 357       for (i = 0; i < nz; i++) {
 358         zstd = zstd + pow(*(z+i) - (zsum/nz),2);
 359       }
 360       zmean = Z2dBZ(zsum/nz);
 361       zstd = sqrt(zstd/(nz-1));
 362       zstd = Z2dBZ(zstd);
 363     }
 364
 365     if ((nv < NMIN) || (nz < NMIN)) {
 366       RaveField_setValue(ff_field, 0, yindex, -9999.0);
 367       RaveField_setValue(ff_dev_field, 0, yindex, -9999.0);
 368       RaveField_setValue(dd_field, 0, yindex, -9999.0);
 369       RaveField_setValue(dbzh_field, 0, yindex, -9999.0);
 370       RaveField_setValue(dbzh_dev_field, 0, yindex, -9999.0);
 371       RaveField_setValue(nz_field, 0, yindex, 0);
 372       RaveField_setValue(nv_field, 0, yindex, 0);
 373     } else {
 374       RaveField_setValue(ff_field, 0, yindex, vvel);
 375       RaveField_setValue(ff_dev_field, 0, yindex, vstd);
 376       RaveField_setValue(dd_field, 0, yindex, vdir);
 377       RaveField_setValue(dbzh_field, 0, yindex, zmean);
 378       RaveField_setValue(dbzh_dev_field, 0, yindex, zstd);
 379       RaveField_setValue(nz_field, 0, yindex, nz);
 380       RaveField_setValue(nv_field, 0, yindex, nv);
 381     }
 382
 383 /*    if ((nv < NMIN) || (nz < NMIN))
 384       printf("%6d %6d %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f \n", iz+self->dz/2, 0, -9999., -9999., -9999., -9999., -9999., -9999.);
 385     else
 386       printf("%6d %6d %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f \n", iz+self->dz/2, nv, vvel, vstd, vdir, -9999., zmean, zstd);*/
 387
 388     RAVE_FREE(A);
 389     RAVE_FREE(b);
 390     RAVE_FREE(v);
 391     RAVE_FREE(z);
 392     RAVE_FREE(az);
 393
 394     yindex++; /* Next interval*/
 395   }
 396
 397   result = RAVE_OBJECT_NEW(&VerticalProfile_TYPE);
 398   if (result != NULL) {
 399     if (!VerticalProfile_setFF(result, ff_field) ||
 400         !VerticalProfile_setFFDev(result, ff_dev_field) ||
 401         !VerticalProfile_setDD(result, dd_field) ||
 402         !VerticalProfile_setDBZ(result, dbzh_field) ||
 403         !VerticalProfile_setDBZDev(result, dbzh_dev_field)) {
 404       RAVE_ERROR0("Failed to set vertical profile fields");
 405       RAVE_OBJECT_RELEASE(result);
 406     }
 407   }
 408   VerticalProfile_setLongitude(result, PolarVolume_getLongitude(inobj));
 409   VerticalProfile_setLatitude(result, PolarVolume_getLatitude(inobj));
 410   VerticalProfile_setHeight(result, PolarVolume_getHeight(inobj));
 411   VerticalProfile_setSource(result, PolarVolume_getSource(inobj));
 412   VerticalProfile_setInterval(result, self->dz);
 413   VerticalProfile_setLevels(result, ysize);
 414   VerticalProfile_setMinheight(result, self->dz / 2);
 415   VerticalProfile_setMaxheight(result, self->hmax);
 416   VerticalProfile_setDate(result, PolarVolume_getDate(inobj));
 417   VerticalProfile_setTime(result, PolarVolume_getTime(inobj));
 418
 419   WrwpInternal_findAndAddAttribute(result, inobj, "how/highprf");
 420   WrwpInternal_findAndAddAttribute(result, inobj, "how/lowprf");
 421   WrwpInternal_findAndAddAttribute(result, inobj, "how/pulsewidth");
 422   WrwpInternal_findAndAddAttribute(result, inobj, "how/wavelength");
 423
 424   WrwpInternal_findAndAddAttribute(result, inobj, "how/RXbandwidth");
 425   WrwpInternal_findAndAddAttribute(result, inobj, "how/RXloss");
 426   WrwpInternal_findAndAddAttribute(result, inobj, "how/TXloss");
 427   WrwpInternal_findAndAddAttribute(result, inobj, "how/antgain");
 428   WrwpInternal_findAndAddAttribute(result, inobj, "how/azmethod");
 429   WrwpInternal_findAndAddAttribute(result, inobj, "how/binmethod");
 430   WrwpInternal_findAndAddAttribute(result, inobj, "how/malfunc");
 431   WrwpInternal_findAndAddAttribute(result, inobj, "how/nomTXpower");
 432   WrwpInternal_findAndAddAttribute(result, inobj, "how/radar_msg");
 433   WrwpInternal_findAndAddAttribute(result, inobj, "how/radconstH");
 434   WrwpInternal_findAndAddAttribute(result, inobj, "how/radomeloss");
 435   WrwpInternal_findAndAddAttribute(result, inobj, "how/rpm");
 436   WrwpInternal_findAndAddAttribute(result, inobj, "how/software");
 437   WrwpInternal_findAndAddAttribute(result, inobj, "how/system");
 438
 439   WrwpInternal_addIntAttribute(result, "how/minrange", Wrwp_getDMIN(self));
 440   WrwpInternal_addIntAttribute(result, "how/maxrange", Wrwp_getDMAX(self));
 441
 442 done:
 443   RAVE_OBJECT_RELEASE(polnav);
 444   RAVE_OBJECT_RELEASE(ff_field);
 445   RAVE_OBJECT_RELEASE(ff_dev_field);
 446   RAVE_OBJECT_RELEASE(dd_field);
 447   RAVE_OBJECT_RELEASE(dbzh_field);
 448   RAVE_OBJECT_RELEASE(dbzh_dev_field);
 449   RAVE_OBJECT_RELEASE(nz_field);
 450   RAVE_OBJECT_RELEASE(nv_field);
 451
 452   return result;
 453 }
 454
 455 /*@} End of Interface functions */
 456
 457 RaveCoreObjectType Wrwp_TYPE = {
 458     "Wrwp",
 459     sizeof(Wrwp_t),
 460     Wrwp_constructor,
 461     Wrwp_destructor
 462 };
 463