Cleaned up code
authorAnders Henja <anders@henjab.se>
Wed, 18 Sep 2013 12:21:18 +0000 (14:21 +0200)
committerAnders Henja <anders@henjab.se>
Wed, 18 Sep 2013 12:21:18 +0000 (14:21 +0200)
lib/wrwp.c
pywrwp/pywrwp.c
test/pytest/WrwpTest.py

index 4eee9e5..af3d661 100644 (file)
@@ -144,154 +144,175 @@ VerticalProfile_t* Wrwp_generate(Wrwp_t* self, PolarVolume_t* inobj) {
        PolarScanParam_t* vrad = NULL;
        PolarScanParam_t* dbz = NULL;
        PolarNavigator_t* polnav = NULL;
-       int nrhs = NRHS, lda = LDA, ldb = LDB, info;
+       int nrhs = NRHS, lda = LDA, ldb = LDB;
        int nscans = 0, nbins = 0, nrays, nv, nz, i, iz, is, ib, ir;
        double rscale, elangle, gain, offset, nodata, undetect, val;
        double d, h;
        double alpha, beta, gamma, vvel, vdir, vstd, zsum, zmean, zstd;
 
+       RAVE_ASSERT((self != NULL), "self == NULL");
+
        result = RAVE_OBJECT_NEW (&VerticalProfile_TYPE);
+
        polnav = RAVE_OBJECT_NEW (&PolarNavigator_TYPE);
        nscans = PolarVolume_getNumberOfScans (inobj);
 
        // loop over atmospheric layers
-       for (iz=0; iz<self->hmax; iz+=self->dz) {
-
-               /* allocate memory and initialize with zeros */
-               double *A = RAVE_CALLOC ((size_t)(NOR*NOC), sizeof (double));
-               double *b = RAVE_CALLOC ((size_t)(NOR), sizeof (double));
-               double *v = RAVE_CALLOC ((size_t)(NOR), sizeof (double));
-               double *z = RAVE_CALLOC ((size_t)(NOR), sizeof (double));
-               double *az = RAVE_CALLOC ((size_t)(NOR), sizeof (double));
-
-               vdir = -9999.0; /*NAN;*/
-               vvel = -9999.0; /*NAN;*/
-               vstd = 0.0;
-               zsum = 0.0;
-               zmean = -9999.0; /*NAN;*/
-               zstd = 0.0;
-               nv = 0;
-               nz = 0;
-
-               // loop over scans
-               for (is=0; is<nscans; is++) {
-                       scan = PolarVolume_getScan (inobj, is);
-                       nbins = PolarScan_getNbins (scan);
-                       nrays = PolarScan_getNrays (scan);
-                       rscale = PolarScan_getRscale (scan);
-                       elangle = PolarScan_getElangle (scan);
-
-                       // radial wind scans
-                       if (PolarScan_hasParameter (scan, "VRAD")) {
-                               vrad = PolarScan_getParameter (scan, "VRAD");
-                               gain = PolarScanParam_getGain (vrad);
-                               offset = PolarScanParam_getOffset (vrad);
-                               nodata = PolarScanParam_getNodata (vrad);
-                               undetect = PolarScanParam_getUndetect (vrad);
-
-                               for (ir=0; ir<nrays; ir++) {
-                                       for (ib=0; ib<nbins; ib++) {
-                                               PolarNavigator_reToDh (polnav, (ib+0.5)*rscale, elangle, &d, &h);
-                                               PolarScanParam_getValue (vrad, ib, ir, &val);
-                                               if ((h>=iz) && (h<iz+self->dz) && (d>=self->dmin) && (d<=self->dmax) && (elangle*RAD2DEG>=self->emin) &&
-                                                       (val!=nodata) && (val!=undetect) && (abs(offset+gain*val)>=self->vmin)) {
-                                                       *(v+nv) = offset+gain*val;
-                                                       *(az+nv) = 360./nrays*ir*DEG2RAD;
-                                                       *(A+nv*NOC) = sin(*(az+nv));
-                                                       *(A+nv*NOC+1) = cos(*(az+nv));
-                                                       *(A+nv*NOC+2) = 1;
-                                                       *(b+nv) = *(v+nv);
-                                                       nv = nv+1;
-                                               }
-                                       }
-                               }
-                       }
-
-                       // reflectivity scans
-                       if (PolarScan_hasParameter (scan, "DBZH")) {
-                               dbz = PolarScan_getParameter (scan, "DBZH");
-                               gain = PolarScanParam_getGain (dbz);
-                               offset = PolarScanParam_getOffset (dbz);
-                               nodata = PolarScanParam_getNodata (dbz);
-                               undetect = PolarScanParam_getUndetect (dbz);
-
-                               for (ir=0; ir<nrays; ir++) {
-                                       for (ib=0; ib<nbins; ib++) {
-                                               PolarNavigator_reToDh (polnav, (ib+0.5)*rscale, elangle, &d, &h);
-                                               PolarScanParam_getValue (dbz, ib, ir, &val);
-                                               if ((h>=iz) && (h<iz+self->dz) && (d>=self->dmin) && (d<=self->dmax) && (elangle*RAD2DEG>=self->emin) &&
-                                                       (val!=nodata) && (val!=undetect)) {
-                                                       *(z+nz) = dBZ2Z (offset+gain*val);
-                                                       zsum = zsum + *(z+nz);
-                                                       nz = nz+1;
-                                               }
-                                       }
-                               }
-                       }
-               }
-
-               // radial wind calculations
-               if (nv>3) {
-
-               //***************************************************************
-               // fitting: y = gamma+alpha*sin(x+beta)                         *
-               // alpha -> amplitude                                           *
-               // beta -> phase shift                                          *
-               // gamma -> consider an y-shift due to the terminal velocity of *
-               //          falling rain drops                                  *
-               //***************************************************************
-
-                       /* QR decomposition */
-                       info = LAPACKE_dgels (LAPACK_ROW_MAJOR, 'N', NOR, NOC, nrhs, A, lda, b, ldb);
-
-                       /* parameter of the wind model */
-                       alpha = sqrt (pow (*(b),2) + pow (*(b+1),2));
-                       beta = atan2 (*(b+1), *b);
-                       gamma = *(b+2);
-
-                       /* wind velocity */
-                       vvel = alpha;
-
-               /* wind direction */
-               vdir = 0;
-               if (alpha < 0) vdir = (PI/2-beta) * RAD2DEG;
-               else if (alpha > 0) vdir = (3*PI/2-beta) * RAD2DEG;
-               if (vdir < 0) vdir = vdir + 360;
-               else if (vdir > 360) vdir = vdir - 360;
-
-               /* standard deviation of the wind velocity */
-               for (i=0; i<nv; i++) {
-                       vstd = vstd + pow (*(v+i) - (gamma+alpha*sin (*(az+i)+beta)),2);
-               }
-               vstd = sqrt (vstd/(nv-1));
-               }
-
-               // reflectivity calculations
-               if (nz>1) {
-
-                       /* standard deviation of reflectivity */
-                       for (i=0; i<nz; i++) {
-                               zstd = zstd + pow (*(z+i) - (zsum/nz),2);
-                       }
-                       zmean = Z2dBZ (zsum/nz);
-                       zstd = sqrt (zstd/(nz-1));
-                       zstd = Z2dBZ (zstd);
-               }
-
-               if ((nv<NMIN) || (nz<NMIN))
-                       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.);
-               else
-                       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);
-
-               RAVE_FREE (A);
-               RAVE_FREE (b);
-               RAVE_FREE (v);
-               RAVE_FREE (z);
-               RAVE_FREE (az);
-       }
+  for (iz = 0; iz < self->hmax; iz += self->dz) {
+    /* allocate memory and initialize with zeros */
+    double *A = RAVE_CALLOC((size_t)(NOR*NOC), sizeof (double));
+    double *b = RAVE_CALLOC((size_t)(NOR), sizeof (double));
+    double *v = RAVE_CALLOC((size_t)(NOR), sizeof (double));
+    double *z = RAVE_CALLOC((size_t)(NOR), sizeof (double));
+    double *az = RAVE_CALLOC((size_t)(NOR), sizeof (double));
+
+    vdir = -9999.0;
+    vvel = -9999.0;
+    vstd = 0.0;
+    zsum = 0.0;
+    zmean = -9999.0;
+    zstd = 0.0;
+    nv = 0;
+    nz = 0;
+
+    // loop over scans
+    for (is = 0; is < nscans; is++) {
+      scan = PolarVolume_getScan(inobj, is);
+      nbins = PolarScan_getNbins(scan);
+      nrays = PolarScan_getNrays(scan);
+      rscale = PolarScan_getRscale(scan);
+      elangle = PolarScan_getElangle(scan);
+
+      // radial wind scans
+      if (PolarScan_hasParameter(scan, "VRAD")) {
+        vrad = PolarScan_getParameter(scan, "VRAD");
+        gain = PolarScanParam_getGain(vrad);
+        offset = PolarScanParam_getOffset(vrad);
+        nodata = PolarScanParam_getNodata(vrad);
+        undetect = PolarScanParam_getUndetect(vrad);
+
+        for (ir = 0; ir < nrays; ir++) {
+          for (ib = 0; ib < nbins; ib++) {
+            PolarNavigator_reToDh(polnav, (ib+0.5)*rscale, elangle, &d, &h);
+            PolarScanParam_getValue(vrad, ib, ir, &val);
+
+            if ((h >= iz) &&
+                (h < iz + self->dz) &&
+                (d >= self->dmin) &&
+                (d <= self->dmax) &&
+                (elangle * RAD2DEG >= self->emin) &&
+                (val != nodata) &&
+                (val != undetect) &&
+                (abs(offset + gain * val) >= self->vmin)) {
+              *(v+nv) = offset+gain*val;
+              *(az+nv) = 360./nrays*ir*DEG2RAD;
+              *(A+nv*NOC) = sin(*(az+nv));
+              *(A+nv*NOC+1) = cos(*(az+nv));
+              *(A+nv*NOC+2) = 1;
+              *(b+nv) = *(v+nv);
+              nv = nv+1;
+            }
+          }
+        }
+        RAVE_OBJECT_RELEASE(vrad);
+      }
+
+      // reflectivity scans
+      if (PolarScan_hasParameter(scan, "DBZH")) {
+        dbz = PolarScan_getParameter(scan, "DBZH");
+        gain = PolarScanParam_getGain(dbz);
+        offset = PolarScanParam_getOffset(dbz);
+        nodata = PolarScanParam_getNodata(dbz);
+        undetect = PolarScanParam_getUndetect(dbz);
+
+        for (ir = 0; ir < nrays; ir++) {
+          for (ib = 0; ib < nbins; ib++) {
+            PolarNavigator_reToDh (polnav, (ib+0.5)*rscale, elangle, &d, &h);
+            PolarScanParam_getValue (dbz, ib, ir, &val);
+            if ((h >= iz) &&
+                (h < iz+self->dz) &&
+                (d >= self->dmin) &&
+                (d <= self->dmax) &&
+                (elangle*RAD2DEG >= self->emin) &&
+                (val != nodata) &&
+                (val != undetect)) {
+              *(z+nz) = dBZ2Z(offset+gain*val);
+              zsum = zsum + *(z+nz);
+              nz = nz+1;
+            }
+          }
+        }
+        RAVE_OBJECT_RELEASE(dbz);
+      }
+      RAVE_OBJECT_RELEASE(scan);
+    }
+
+    // radial wind calculations
+    if (nv>3) {
+      //***************************************************************
+      // fitting: y = gamma+alpha*sin(x+beta)                         *
+      // alpha -> amplitude                                           *
+      // beta -> phase shift                                          *
+      // gamma -> consider an y-shift due to the terminal velocity of *
+      //          falling rain drops                                  *
+      //***************************************************************
+
+      /* QR decomposition */
+      /*info = */LAPACKE_dgels(LAPACK_ROW_MAJOR, 'N', NOR, NOC, nrhs, A, lda, b, ldb);
+
+      /* parameter of the wind model */
+      alpha = sqrt(pow(*(b),2) + pow(*(b+1),2));
+      beta = atan2(*(b+1), *b);
+      gamma = *(b+2);
+
+      /* wind velocity */
+      vvel = alpha;
+
+      /* wind direction */
+      vdir = 0;
+      if (alpha < 0) {
+        vdir = (PI/2-beta) * RAD2DEG;
+      } else if (alpha > 0) {
+        vdir = (3*PI/2-beta) * RAD2DEG;
+      }
+      if (vdir < 0) {
+        vdir = vdir + 360;
+      } else if (vdir > 360) {
+        vdir = vdir - 360;
+      }
+
+      /* standard deviation of the wind velocity */
+      for (i=0; i<nv; i++) {
+        vstd = vstd + pow (*(v+i) - (gamma+alpha*sin (*(az+i)+beta)),2);
+      }
+      vstd = sqrt (vstd/(nv-1));
+    }
+
+    // reflectivity calculations
+    if (nz > 1) {
+      /* standard deviation of reflectivity */
+      for (i = 0; i < nz; i++) {
+        zstd = zstd + pow(*(z+i) - (zsum/nz),2);
+      }
+      zmean = Z2dBZ(zsum/nz);
+      zstd = sqrt(zstd/(nz-1));
+      zstd = Z2dBZ(zstd);
+    }
+
+    if ((nv < NMIN) || (nz < NMIN))
+      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.);
+    else
+      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);
+
+    RAVE_FREE(A);
+    RAVE_FREE(b);
+    RAVE_FREE(v);
+    RAVE_FREE(z);
+    RAVE_FREE(az);
+  }
 
   RAVE_OBJECT_RELEASE(polnav);
-       return result;
+  return result;
 }
 
 /*@} End of Interface functions */
index aa3be93..4e401fb 100644 (file)
@@ -150,16 +150,28 @@ static PyObject* _pywrwp_new(PyObject* self, PyObject* args)
 static PyObject* _pywrwp_generate(PyWrwp* self, PyObject* args)
 {
        PyObject* obj = NULL;
-       //PyPolarVolume* pyobj = NULL;
+       PyVerticalProfile* pyvp = NULL;
+       VerticalProfile_t* vp = NULL;
+
        if(!PyArg_ParseTuple(args, "O",&obj)) {
                return NULL;
        }
+
        if (!PyPolarVolume_Check(obj)) {
          raiseException_returnNULL(PyExc_AttributeError, "In argument must be a polar volume");
        }
-       //pyobj = (PyPolarVolume*)obj;
 
-       return NULL;
+       vp = Wrwp_generate(self->wrwp, ((PyPolarVolume*)obj)->pvol);
+
+       if (vp == NULL) {
+         raiseException_gotoTag(done, PyExc_RuntimeError, "Failed to generate vertical profile");
+       }
+
+       pyvp = PyVerticalProfile_New(vp);
+
+done:
+  RAVE_OBJECT_RELEASE(vp);
+       return (PyObject*)pyvp;
 }
 
 /**
@@ -168,14 +180,6 @@ static PyObject* _pywrwp_generate(PyWrwp* self, PyObject* args)
 static struct PyMethodDef _pywrwp_methods[] =
 {
   {"generate", (PyCFunction)_pywrwp_generate, 1},
-/*
-  {"method", NULL},
-  {"ppi", (PyCFunction) _pytransform_ppi, 1},
-  {"cappi", (PyCFunction) _pytransform_cappi, 1},
-  {"pcappi", (PyCFunction) _pytransform_pcappi, 1},
-  {"ctoscan", (PyCFunction) _pytransform_ctoscan, 1},
-  {"ctop", (PyCFunction) _pytransform_ctop, 1},
-  {"fillGap", (PyCFunction) _pytransform_fillGap, 1},*/
   {NULL, NULL } /* sentinel */
 };
 
@@ -326,6 +330,7 @@ init_wrwp(void)
 
   import_array();
   import_pypolarvolume();
+  import_pyverticalprofile();
   PYRAVE_DEBUG_INITIALIZE;
 }
 /*@} End of Module setup */
index 1917bb4..919287e 100644 (file)
@@ -28,8 +28,11 @@ import unittest
 import string
 import _wrwp
 import _helpers
+import _raveio, _rave
 
 class WrwpTest(unittest.TestCase):
+  FIXTURE = "fixtures/pvol_seang_20090501T120000Z.h5"
+  
   def setUp(self):
     _helpers.triggerMemoryStatus()
 
@@ -103,6 +106,11 @@ class WrwpTest(unittest.TestCase):
     self.assertAlmostEquals(3.5, obj.vmin, 4)
     obj.vmin = 4
     self.assertAlmostEquals(4.0, obj.vmin, 4)
+  
+  def test_generate(self):
+    pvol = _raveio.open(self.FIXTURE).object
+    generator = _wrwp.new()
+    vp = generator.generate(pvol)
     
 if __name__ == "__main__":
   unittest.main()
\ No newline at end of file