Cleaned up code
[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
153         RAVE_ASSERT((self != NULL), "self == NULL");
154
155         result = RAVE_OBJECT_NEW (&VerticalProfile_TYPE);
156
157         polnav = RAVE_OBJECT_NEW (&PolarNavigator_TYPE);
158         nscans = PolarVolume_getNumberOfScans (inobj);
159
160         // loop over atmospheric layers
161   for (iz = 0; iz < self->hmax; iz += self->dz) {
162     /* allocate memory and initialize with zeros */
163     double *A = RAVE_CALLOC((size_t)(NOR*NOC), sizeof (double));
164     double *b = RAVE_CALLOC((size_t)(NOR), sizeof (double));
165     double *v = RAVE_CALLOC((size_t)(NOR), sizeof (double));
166     double *z = RAVE_CALLOC((size_t)(NOR), sizeof (double));
167     double *az = RAVE_CALLOC((size_t)(NOR), sizeof (double));
168
169     vdir = -9999.0;
170     vvel = -9999.0;
171     vstd = 0.0;
172     zsum = 0.0;
173     zmean = -9999.0;
174     zstd = 0.0;
175     nv = 0;
176     nz = 0;
177
178     // loop over scans
179     for (is = 0; is < nscans; is++) {
180       scan = PolarVolume_getScan(inobj, is);
181       nbins = PolarScan_getNbins(scan);
182       nrays = PolarScan_getNrays(scan);
183       rscale = PolarScan_getRscale(scan);
184       elangle = PolarScan_getElangle(scan);
185
186       // radial wind scans
187       if (PolarScan_hasParameter(scan, "VRAD")) {
188         vrad = PolarScan_getParameter(scan, "VRAD");
189         gain = PolarScanParam_getGain(vrad);
190         offset = PolarScanParam_getOffset(vrad);
191         nodata = PolarScanParam_getNodata(vrad);
192         undetect = PolarScanParam_getUndetect(vrad);
193
194         for (ir = 0; ir < nrays; ir++) {
195           for (ib = 0; ib < nbins; ib++) {
196             PolarNavigator_reToDh(polnav, (ib+0.5)*rscale, elangle, &d, &h);
197             PolarScanParam_getValue(vrad, ib, ir, &val);
198
199             if ((h >= iz) &&
200                 (h < iz + self->dz) &&
201                 (d >= self->dmin) &&
202                 (d <= self->dmax) &&
203                 (elangle * RAD2DEG >= self->emin) &&
204                 (val != nodata) &&
205                 (val != undetect) &&
206                 (abs(offset + gain * val) >= self->vmin)) {
207               *(v+nv) = offset+gain*val;
208               *(az+nv) = 360./nrays*ir*DEG2RAD;
209               *(A+nv*NOC) = sin(*(az+nv));
210               *(A+nv*NOC+1) = cos(*(az+nv));
211               *(A+nv*NOC+2) = 1;
212               *(b+nv) = *(v+nv);
213               nv = nv+1;
214             }
215           }
216         }
217         RAVE_OBJECT_RELEASE(vrad);
218       }
219
220       // reflectivity scans
221       if (PolarScan_hasParameter(scan, "DBZH")) {
222         dbz = PolarScan_getParameter(scan, "DBZH");
223         gain = PolarScanParam_getGain(dbz);
224         offset = PolarScanParam_getOffset(dbz);
225         nodata = PolarScanParam_getNodata(dbz);
226         undetect = PolarScanParam_getUndetect(dbz);
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 (dbz, ib, ir, &val);
232             if ((h >= iz) &&
233                 (h < iz+self->dz) &&
234                 (d >= self->dmin) &&
235                 (d <= self->dmax) &&
236                 (elangle*RAD2DEG >= self->emin) &&
237                 (val != nodata) &&
238                 (val != undetect)) {
239               *(z+nz) = dBZ2Z(offset+gain*val);
240               zsum = zsum + *(z+nz);
241               nz = nz+1;
242             }
243           }
244         }
245         RAVE_OBJECT_RELEASE(dbz);
246       }
247       RAVE_OBJECT_RELEASE(scan);
248     }
249
250     // radial wind calculations
251     if (nv>3) {
252       //***************************************************************
253       // fitting: y = gamma+alpha*sin(x+beta)                         *
254       // alpha -> amplitude                                           *
255       // beta -> phase shift                                          *
256       // gamma -> consider an y-shift due to the terminal velocity of *
257       //          falling rain drops                                  *
258       //***************************************************************
259
260       /* QR decomposition */
261       /*info = */LAPACKE_dgels(LAPACK_ROW_MAJOR, 'N', NOR, NOC, nrhs, A, lda, b, ldb);
262
263       /* parameter of the wind model */
264       alpha = sqrt(pow(*(b),2) + pow(*(b+1),2));
265       beta = atan2(*(b+1), *b);
266       gamma = *(b+2);
267
268       /* wind velocity */
269       vvel = alpha;
270
271       /* wind direction */
272       vdir = 0;
273       if (alpha < 0) {
274         vdir = (PI/2-beta) * RAD2DEG;
275       } else if (alpha > 0) {
276         vdir = (3*PI/2-beta) * RAD2DEG;
277       }
278       if (vdir < 0) {
279         vdir = vdir + 360;
280       } else if (vdir > 360) {
281         vdir = vdir - 360;
282       }
283
284       /* standard deviation of the wind velocity */
285       for (i=0; i<nv; i++) {
286         vstd = vstd + pow (*(v+i) - (gamma+alpha*sin (*(az+i)+beta)),2);
287       }
288       vstd = sqrt (vstd/(nv-1));
289     }
290
291     // reflectivity calculations
292     if (nz > 1) {
293       /* standard deviation of reflectivity */
294       for (i = 0; i < nz; i++) {
295         zstd = zstd + pow(*(z+i) - (zsum/nz),2);
296       }
297       zmean = Z2dBZ(zsum/nz);
298       zstd = sqrt(zstd/(nz-1));
299       zstd = Z2dBZ(zstd);
300     }
301
302     if ((nv < NMIN) || (nz < NMIN))
303       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.);
304     else
305       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);
306
307     RAVE_FREE(A);
308     RAVE_FREE(b);
309     RAVE_FREE(v);
310     RAVE_FREE(z);
311     RAVE_FREE(az);
312   }
313
314   RAVE_OBJECT_RELEASE(polnav);
315   return result;
316 }
317
318 /*@} End of Interface functions */
319
320 RaveCoreObjectType Wrwp_TYPE = {
321     "Wrwp",
322     sizeof(Wrwp_t),
323     Wrwp_constructor,
324     Wrwp_destructor
325 };
326