Ticket 504: How attributes should be copied from volume to profile
[baltrad-wrwp.git] / 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