7f831f426e9e029455a4bf1d52a3948ad45e530b
[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 /*@} End of Private functions */
67
68 /*@{ Interface functions */
69 void Wrwp_setDZ(Wrwp_t* self, int dz)
70 {
71   RAVE_ASSERT((self != NULL), "self == NULL");
72   self->dz = dz;
73 }
74
75 int Wrwp_getDZ(Wrwp_t* self)
76 {
77   RAVE_ASSERT((self != NULL), "self == NULL");
78   return self->dz;
79 }
80
81 void Wrwp_setHMAX(Wrwp_t* self, int hmax)
82 {
83   RAVE_ASSERT((self != NULL), "self == NULL");
84   self->hmax = hmax;
85 }
86
87 int Wrwp_getHMAX(Wrwp_t* self)
88 {
89   RAVE_ASSERT((self != NULL), "self == NULL");
90   return self->hmax;
91 }
92
93 void Wrwp_setDMIN(Wrwp_t* self, int dmin)
94 {
95   RAVE_ASSERT((self != NULL), "self == NULL");
96   self->dmin = dmin;
97 }
98
99 int Wrwp_getDMIN(Wrwp_t* self)
100 {
101   RAVE_ASSERT((self != NULL), "self == NULL");
102   return self->dmin;
103 }
104
105 void Wrwp_setDMAX(Wrwp_t* self, int dmax)
106 {
107   RAVE_ASSERT((self != NULL), "self == NULL");
108   self->dmax = dmax;
109 }
110
111 int Wrwp_getDMAX(Wrwp_t* self)
112 {
113   RAVE_ASSERT((self != NULL), "self == NULL");
114   return self->dmax;
115 }
116
117 void Wrwp_setEMIN(Wrwp_t* self, double emin)
118 {
119   RAVE_ASSERT((self != NULL), "self == NULL");
120   self->emin = emin;
121 }
122
123 double Wrwp_getEMIN(Wrwp_t* self)
124 {
125   RAVE_ASSERT((self != NULL), "self == NULL");
126   return self->emin;
127 }
128
129 void Wrwp_setVMIN(Wrwp_t* self, double vmin)
130 {
131   RAVE_ASSERT((self != NULL), "self == NULL");
132   self->vmin = vmin;
133 }
134
135 double Wrwp_getVMIN(Wrwp_t* self)
136 {
137   RAVE_ASSERT((self != NULL), "self == NULL");
138   return self->vmin;
139 }
140
141 VerticalProfile_t* Wrwp_generate(Wrwp_t* self, PolarVolume_t* inobj) {
142   VerticalProfile_t* result = NULL;
143         PolarScan_t* scan = NULL;
144         PolarScanParam_t* vrad = NULL;
145         PolarScanParam_t* dbz = NULL;
146         PolarNavigator_t* polnav = NULL;
147         int nrhs = NRHS, lda = LDA, ldb = LDB;
148         int nscans = 0, nbins = 0, nrays, nv, nz, i, iz, is, ib, ir;
149         double rscale, elangle, gain, offset, nodata, undetect, val;
150         double d, h;
151         double alpha, beta, gamma, vvel, vdir, vstd, zsum, zmean, zstd;
152         int ysize = 0, yindex = 0;
153
154         RaveField_t *ff_field = NULL, *ff_dev_field = NULL, *dd_field = NULL;
155         RaveField_t *dbzh_field = NULL, *dbzh_dev_field = NULL, *nz_field = NULL, *nv_field = NULL;
156
157         RAVE_ASSERT((self != NULL), "self == NULL");
158
159         ff_field = RAVE_OBJECT_NEW(&RaveField_TYPE);
160   ff_dev_field = RAVE_OBJECT_NEW(&RaveField_TYPE);
161   dd_field = RAVE_OBJECT_NEW(&RaveField_TYPE);
162   dbzh_field = RAVE_OBJECT_NEW(&RaveField_TYPE);
163   dbzh_dev_field = RAVE_OBJECT_NEW(&RaveField_TYPE);
164   nz_field = RAVE_OBJECT_NEW(&RaveField_TYPE);
165   nv_field = RAVE_OBJECT_NEW(&RaveField_TYPE);
166
167   if (ff_field == NULL || ff_dev_field == NULL || dd_field == NULL ||
168       dbzh_field == NULL || dbzh_dev_field == NULL || nz_field == NULL || nv_field == NULL) {
169     RAVE_ERROR0("Failed to allocate memory for the resulting vp fields");
170     goto done;
171   }
172   ysize = self->hmax / self->dz;
173   if (!RaveField_createData(ff_field, 1, ysize, RaveDataType_DOUBLE) ||
174       !RaveField_createData(ff_dev_field, 1, ysize, RaveDataType_DOUBLE) ||
175       !RaveField_createData(dd_field, 1, ysize, RaveDataType_DOUBLE) ||
176       !RaveField_createData(dbzh_field, 1, ysize, RaveDataType_DOUBLE) ||
177       !RaveField_createData(dbzh_dev_field, 1, ysize, RaveDataType_DOUBLE) ||
178       !RaveField_createData(nz_field, 1, ysize, RaveDataType_INT) ||
179       !RaveField_createData(nv_field, 1, ysize, RaveDataType_INT)) {
180     RAVE_ERROR0("Failed to allocate arrays for the resulting vp fields");
181     goto done;
182   }
183
184         polnav = RAVE_OBJECT_NEW (&PolarNavigator_TYPE);
185         PolarNavigator_setLat0(polnav, PolarVolume_getLatitude(inobj));
186         PolarNavigator_setLon0(polnav, PolarVolume_getLongitude(inobj));
187         PolarNavigator_setAlt0(polnav, PolarVolume_getHeight(inobj));
188
189         nscans = PolarVolume_getNumberOfScans (inobj);
190
191         // We use yindex for filling in the arrays even though we loop to hmax...
192         yindex = 0;
193
194         // loop over atmospheric layers
195   for (iz = 0; iz < self->hmax; iz += self->dz) {
196     /* allocate memory and initialize with zeros */
197     double *A = RAVE_CALLOC((size_t)(NOR*NOC), sizeof (double));
198     double *b = RAVE_CALLOC((size_t)(NOR), sizeof (double));
199     double *v = RAVE_CALLOC((size_t)(NOR), sizeof (double));
200     double *z = RAVE_CALLOC((size_t)(NOR), sizeof (double));
201     double *az = RAVE_CALLOC((size_t)(NOR), sizeof (double));
202
203     vdir = -9999.0;
204     vvel = -9999.0;
205     vstd = 0.0;
206     zsum = 0.0;
207     zmean = -9999.0;
208     zstd = 0.0;
209     nv = 0;
210     nz = 0;
211
212     // loop over scans
213     for (is = 0; is < nscans; is++) {
214       scan = PolarVolume_getScan(inobj, is);
215       nbins = PolarScan_getNbins(scan);
216       nrays = PolarScan_getNrays(scan);
217       rscale = PolarScan_getRscale(scan);
218       elangle = PolarScan_getElangle(scan);
219
220       // radial wind scans
221       if (PolarScan_hasParameter(scan, "VRAD")) {
222         vrad = PolarScan_getParameter(scan, "VRAD");
223         gain = PolarScanParam_getGain(vrad);
224         offset = PolarScanParam_getOffset(vrad);
225         nodata = PolarScanParam_getNodata(vrad);
226         undetect = PolarScanParam_getUndetect(vrad);
227
228         for (ir = 0; ir < nrays; ir++) {
229           for (ib = 0; ib < nbins; ib++) {
230             PolarNavigator_reToDh(polnav, (ib+0.5)*rscale, elangle, &d, &h);
231             PolarScanParam_getValue(vrad, ib, ir, &val);
232
233             if ((h >= iz) &&
234                 (h < iz + self->dz) &&
235                 (d >= self->dmin) &&
236                 (d <= self->dmax) &&
237                 (elangle * RAD2DEG >= self->emin) &&
238                 (val != nodata) &&
239                 (val != undetect) &&
240                 (abs(offset + gain * val) >= self->vmin)) {
241               *(v+nv) = offset+gain*val;
242               *(az+nv) = 360./nrays*ir*DEG2RAD;
243               *(A+nv*NOC) = sin(*(az+nv));
244               *(A+nv*NOC+1) = cos(*(az+nv));
245               *(A+nv*NOC+2) = 1;
246               *(b+nv) = *(v+nv);
247               nv = nv+1;
248             }
249           }
250         }
251         RAVE_OBJECT_RELEASE(vrad);
252       }
253
254       // reflectivity scans
255       if (PolarScan_hasParameter(scan, "DBZH")) {
256         dbz = PolarScan_getParameter(scan, "DBZH");
257         gain = PolarScanParam_getGain(dbz);
258         offset = PolarScanParam_getOffset(dbz);
259         nodata = PolarScanParam_getNodata(dbz);
260         undetect = PolarScanParam_getUndetect(dbz);
261
262         for (ir = 0; ir < nrays; ir++) {
263           for (ib = 0; ib < nbins; ib++) {
264             PolarNavigator_reToDh (polnav, (ib+0.5)*rscale, elangle, &d, &h);
265             PolarScanParam_getValue (dbz, ib, ir, &val);
266             if ((h >= iz) &&
267                 (h < iz+self->dz) &&
268                 (d >= self->dmin) &&
269                 (d <= self->dmax) &&
270                 (elangle*RAD2DEG >= self->emin) &&
271                 (val != nodata) &&
272                 (val != undetect)) {
273               *(z+nz) = dBZ2Z(offset+gain*val);
274               zsum = zsum + *(z+nz);
275               nz = nz+1;
276             }
277           }
278         }
279         RAVE_OBJECT_RELEASE(dbz);
280       }
281       RAVE_OBJECT_RELEASE(scan);
282     }
283
284     // radial wind calculations
285     if (nv>3) {
286       //***************************************************************
287       // fitting: y = gamma+alpha*sin(x+beta)                         *
288       // alpha -> amplitude                                           *
289       // beta -> phase shift                                          *
290       // gamma -> consider an y-shift due to the terminal velocity of *
291       //          falling rain drops                                  *
292       //***************************************************************
293
294       /* QR decomposition */
295       /*info = */LAPACKE_dgels(LAPACK_ROW_MAJOR, 'N', NOR, NOC, nrhs, A, lda, b, ldb);
296
297       /* parameter of the wind model */
298       alpha = sqrt(pow(*(b),2) + pow(*(b+1),2));
299       beta = atan2(*(b+1), *b);
300       gamma = *(b+2);
301
302       /* wind velocity */
303       vvel = alpha;
304
305       /* wind direction */
306       vdir = 0;
307       if (alpha < 0) {
308         vdir = (PI/2-beta) * RAD2DEG;
309       } else if (alpha > 0) {
310         vdir = (3*PI/2-beta) * RAD2DEG;
311       }
312       if (vdir < 0) {
313         vdir = vdir + 360;
314       } else if (vdir > 360) {
315         vdir = vdir - 360;
316       }
317
318       /* standard deviation of the wind velocity */
319       for (i=0; i<nv; i++) {
320         vstd = vstd + pow (*(v+i) - (gamma+alpha*sin (*(az+i)+beta)),2);
321       }
322       vstd = sqrt (vstd/(nv-1));
323     }
324
325     // reflectivity calculations
326     if (nz > 1) {
327       /* standard deviation of reflectivity */
328       for (i = 0; i < nz; i++) {
329         zstd = zstd + pow(*(z+i) - (zsum/nz),2);
330       }
331       zmean = Z2dBZ(zsum/nz);
332       zstd = sqrt(zstd/(nz-1));
333       zstd = Z2dBZ(zstd);
334     }
335
336     if ((nv < NMIN) || (nz < NMIN)) {
337       RaveField_setValue(ff_field, 0, yindex, -9999.0);
338       RaveField_setValue(ff_dev_field, 0, yindex, -9999.0);
339       RaveField_setValue(dd_field, 0, yindex, -9999.0);
340       RaveField_setValue(dbzh_field, 0, yindex, -9999.0);
341       RaveField_setValue(dbzh_dev_field, 0, yindex, -9999.0);
342       RaveField_setValue(nz_field, 0, yindex, 0);
343       RaveField_setValue(nv_field, 0, yindex, 0);
344     } else {
345       RaveField_setValue(ff_field, 0, yindex, vvel);
346       RaveField_setValue(ff_dev_field, 0, yindex, vstd);
347       RaveField_setValue(dd_field, 0, yindex, vdir);
348       RaveField_setValue(dbzh_field, 0, yindex, zmean);
349       RaveField_setValue(dbzh_dev_field, 0, yindex, zstd);
350       RaveField_setValue(nz_field, 0, yindex, nz);
351       RaveField_setValue(nv_field, 0, yindex, nv);
352     }
353
354 /*    if ((nv < NMIN) || (nz < NMIN))
355       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.);
356     else
357       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);*/
358
359     RAVE_FREE(A);
360     RAVE_FREE(b);
361     RAVE_FREE(v);
362     RAVE_FREE(z);
363     RAVE_FREE(az);
364
365     yindex++; /* Next interval*/
366   }
367
368   result = RAVE_OBJECT_NEW(&VerticalProfile_TYPE);
369   if (result != NULL) {
370     if (!VerticalProfile_setFF(result, ff_field) ||
371         !VerticalProfile_setFFDev(result, ff_dev_field) ||
372         !VerticalProfile_setDD(result, dd_field) ||
373         !VerticalProfile_setDBZ(result, dbzh_field) ||
374         !VerticalProfile_setDBZDev(result, dbzh_dev_field)) {
375       RAVE_ERROR0("Failed to set vertical profile fields");
376       RAVE_OBJECT_RELEASE(result);
377     }
378   }
379   VerticalProfile_setLongitude(result, PolarVolume_getLongitude(inobj));
380   VerticalProfile_setLatitude(result, PolarVolume_getLatitude(inobj));
381   VerticalProfile_setHeight(result, PolarVolume_getHeight(inobj));
382   VerticalProfile_setSource(result, PolarVolume_getSource(inobj));
383   VerticalProfile_setInterval(result, self->dz);
384   VerticalProfile_setLevels(result, ysize);
385   VerticalProfile_setMinheight(result, self->dz / 2);
386   VerticalProfile_setMaxheight(result, self->hmax);
387   VerticalProfile_setDate(result, PolarVolume_getDate(inobj));
388   VerticalProfile_setTime(result, PolarVolume_getTime(inobj));
389
390 done:
391   RAVE_OBJECT_RELEASE(polnav);
392   RAVE_OBJECT_RELEASE(ff_field);
393   RAVE_OBJECT_RELEASE(ff_dev_field);
394   RAVE_OBJECT_RELEASE(dd_field);
395   RAVE_OBJECT_RELEASE(dbzh_field);
396   RAVE_OBJECT_RELEASE(dbzh_dev_field);
397   RAVE_OBJECT_RELEASE(nz_field);
398   RAVE_OBJECT_RELEASE(nv_field);
399
400   return result;
401 }
402
403 /*@} End of Interface functions */
404
405 RaveCoreObjectType Wrwp_TYPE = {
406     "Wrwp",
407     sizeof(Wrwp_t),
408     Wrwp_constructor,
409     Wrwp_destructor
410 };
411