1 /* --------------------------------------------------------------------
2 Copyright (C) 2011 Swedish Meteorological and Hydrological Institute, SMHI
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.
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.
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 ------------------------------------------------------------------------*/
18 /** Function for deriving weather radar wind and reflectivity profiles
21 * @author Gunther Haase, SMHI
26 #include "rave_debug.h"
27 #include "rave_alloc.h"
31 * Represents one wrwp generator
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] */
43 /*@{ Private functions */
47 static int Wrwp_constructor(RaveCoreObject* obj)
49 Wrwp_t* wrwp = (Wrwp_t*)obj;
62 static void Wrwp_destructor(RaveCoreObject* obj)
66 static int WrwpInternal_findAndAddAttribute(VerticalProfile_t* vp, PolarVolume_t* pvol, const char* name)
68 int nscans = PolarVolume_getNumberOfScans(pvol);
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);
77 RAVE_OBJECT_RELEASE(attr);
79 RAVE_OBJECT_RELEASE(scan);
84 static int WrwpInternal_addIntAttribute(VerticalProfile_t* vp, const char* name, int value)
86 RaveAttribute_t* attr = RaveAttributeHelp_createLong(name, value);
89 result = VerticalProfile_addAttribute(vp, attr);
91 RAVE_OBJECT_RELEASE(attr);
95 /*@} End of Private functions */
97 /*@{ Interface functions */
98 void Wrwp_setDZ(Wrwp_t* self, int dz)
100 RAVE_ASSERT((self != NULL), "self == NULL");
104 int Wrwp_getDZ(Wrwp_t* self)
106 RAVE_ASSERT((self != NULL), "self == NULL");
110 void Wrwp_setHMAX(Wrwp_t* self, int hmax)
112 RAVE_ASSERT((self != NULL), "self == NULL");
116 int Wrwp_getHMAX(Wrwp_t* self)
118 RAVE_ASSERT((self != NULL), "self == NULL");
122 void Wrwp_setDMIN(Wrwp_t* self, int dmin)
124 RAVE_ASSERT((self != NULL), "self == NULL");
128 int Wrwp_getDMIN(Wrwp_t* self)
130 RAVE_ASSERT((self != NULL), "self == NULL");
134 void Wrwp_setDMAX(Wrwp_t* self, int dmax)
136 RAVE_ASSERT((self != NULL), "self == NULL");
140 int Wrwp_getDMAX(Wrwp_t* self)
142 RAVE_ASSERT((self != NULL), "self == NULL");
146 void Wrwp_setEMIN(Wrwp_t* self, double emin)
148 RAVE_ASSERT((self != NULL), "self == NULL");
152 double Wrwp_getEMIN(Wrwp_t* self)
154 RAVE_ASSERT((self != NULL), "self == NULL");
158 void Wrwp_setVMIN(Wrwp_t* self, double vmin)
160 RAVE_ASSERT((self != NULL), "self == NULL");
164 double Wrwp_getVMIN(Wrwp_t* self)
166 RAVE_ASSERT((self != NULL), "self == NULL");
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;
180 double alpha, beta, gamma, vvel, vdir, vstd, zsum, zmean, zstd;
181 int ysize = 0, yindex = 0;
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;
186 RAVE_ASSERT((self != NULL), "self == NULL");
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);
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");
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");
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));
218 nscans = PolarVolume_getNumberOfScans (inobj);
220 // We use yindex for filling in the arrays even though we loop to hmax...
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));
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);
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);
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);
263 (h < iz + self->dz) &&
266 (elangle * RAD2DEG >= self->emin) &&
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));
280 RAVE_OBJECT_RELEASE(vrad);
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);
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);
299 (elangle*RAD2DEG >= self->emin) &&
302 *(z+nz) = dBZ2Z(offset+gain*val);
303 zsum = zsum + *(z+nz);
308 RAVE_OBJECT_RELEASE(dbz);
310 RAVE_OBJECT_RELEASE(scan);
313 // radial wind calculations
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 //***************************************************************
323 /* QR decomposition */
324 /*info = */LAPACKE_dgels(LAPACK_ROW_MAJOR, 'N', NOR, NOC, nrhs, A, lda, b, ldb);
326 /* parameter of the wind model */
327 alpha = sqrt(pow(*(b),2) + pow(*(b+1),2));
328 beta = atan2(*(b+1), *b);
337 vdir = (PI/2-beta) * RAD2DEG;
338 } else if (alpha > 0) {
339 vdir = (3*PI/2-beta) * RAD2DEG;
343 } else if (vdir > 360) {
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);
351 vstd = sqrt (vstd/(nv-1));
354 // reflectivity calculations
356 /* standard deviation of reflectivity */
357 for (i = 0; i < nz; i++) {
358 zstd = zstd + pow(*(z+i) - (zsum/nz),2);
360 zmean = Z2dBZ(zsum/nz);
361 zstd = sqrt(zstd/(nz-1));
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);
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);
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.);
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);*/
394 yindex++; /* Next interval*/
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);
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));
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");
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");
439 WrwpInternal_addIntAttribute(result, "how/minrange", Wrwp_getDMIN(self));
440 WrwpInternal_addIntAttribute(result, "how/maxrange", Wrwp_getDMAX(self));
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);
455 /*@} End of Interface functions */
457 RaveCoreObjectType Wrwp_TYPE = {