First draft of beam blockage hudson-beamb-4-SUCCESS hudson-beamb-5-SUCCESS hudson-beamb-6-SUCCESS hudson-beamb-7-SUCCESS hudson-beamb-8-SUCCESS
authorAnders Henja <anders@baltrad.eu>
Wed, 23 Nov 2011 14:53:37 +0000 (15:53 +0100)
committerAnders Henja <anders@baltrad.eu>
Wed, 23 Nov 2011 14:53:37 +0000 (15:53 +0100)
lib/bbtopography.c
lib/bbtopography.h
lib/beamblockage.c
lib/beamblockage.h
lib/beamblockagemap.c
pybeamb/pybeamblockage.c
pybeamb/pybeamblockagemap.c
test/pytest/PyBeamBlockageMapTest.py
test/pytest/PyBeamBlockageTest.py
test/pytest/fixtures/scan_sevil_20100702T113200Z.h5 [new file with mode: 0644]

index d32862a..55cf1da 100644 (file)
@@ -246,6 +246,31 @@ int BBTopography_setValue(BBTopography_t* self, long col, long row, double value
   return RaveData2D_setValue(self->data, col, row, value);
 }
 
+int BBTopography_getValueAtLonLat(BBTopography_t* self, double lon, double lat, double* v)
+{
+  int result = 0;
+  long ci = 0, ri = 0;
+  double nv = 0.0;
+
+  RAVE_ASSERT((self != NULL), "self == NULL");
+  RAVE_ASSERT((v != NULL), "v == NULL");
+  *v = self->nodata;
+
+  if (self->xdim == 0.0 || self->ydim == 0.0) {
+    RAVE_CRITICAL0("xdim or ydim == 0.0 in topography field");
+    goto done;
+  }
+
+  ci = (lon - self->ulxmap)/self->xdim;
+  ri = (self->ulymap - lat)/self->ydim;
+  if (RaveData2D_getValue(self->data, ci, ri, &nv)) {
+    *v = nv;
+    result = 1;
+  }
+
+done:
+  return result;
+}
 BBTopography_t* BBTopography_concatX(BBTopography_t* self, BBTopography_t* other)
 {
   BBTopography_t *result = NULL;
index c1ef68f..a3ee7ba 100644 (file)
@@ -193,6 +193,19 @@ int BBTopography_getValue(BBTopography_t* self, long col, long row, double* v);
 int BBTopography_setValue(BBTopography_t* self, long col, long row, double value);
 
 /**
+ * Returns the value at the specified lon/lat coordinate. If outside boundaries or if
+ * there is no data at provided coordinate the returned value will be 0 but v will still
+ * be set to nodata.
+ *
+ * @param[in] self - self
+ * @param[in] lon - the longitude
+ * @param[in] lat - the latitude
+ * @param[out] v - the value at specified position
+ * @return the outcome of the operation.
+ */
+int BBTopography_getValueAtLonLat(BBTopography_t* self, double lon, double lat, double* v);
+
+/**
  * Concatenates two topography fields horizontally with each other.
  * The fields and others y-dimension must be the same as well as the data
  * type. It is also necessary that ydim and xdim are the same. All other
index 58d1b6b..0bc2b34 100644 (file)
@@ -130,6 +130,32 @@ done:
   return result;
 }
 
+/**
+ * Generate a running cumulative maximum.
+ * @param[in] p - pointer to array of values
+ * @param[in] nlen - length of vector p
+ */
+static void BeamBlockageInternal_cummax(double* p, long nlen)
+{
+  double t = 0.0;
+  long i = 0;
+
+  RAVE_ASSERT((p != NULL), "p == NULL");
+  if (nlen <= 0) {
+    RAVE_WARNING0("Trying to generate a cumulative max without any data");
+    return;
+  }
+
+  t = p[0];
+  for (i = 1; i < nlen; i++) {
+    if (p[i] < t) {
+      p[i] = t;
+    } else {
+      t = p[i];
+    }
+  }
+}
+
 /*@} End of Private functions */
 
 /*@{ Interface functions */
@@ -145,56 +171,104 @@ const char* BeamBlockage_getTopo30Directory(BeamBlockage_t* self)
   return (const char*)BeamBlockageMap_getTopo30Directory(self->mapper);
 }
 
-RaveField_t* BeamBlockage_getBlockage(BeamBlockage_t* self, PolarScan_t* scan)
+RaveField_t* BeamBlockage_getBlockage(BeamBlockage_t* self, PolarScan_t* scan, double dBlim)
 {
   RaveField_t *field = NULL, *result = NULL;
-  RaveData2D_t *lon = NULL, *lat = NULL;
-  BBTopography_t* topo = NULL;
-  double lon_min = 0.0, lon_max = 0.0, lat_min = 0.0, lat_max = 0.0;
-
-  long ulx = 0, uly = 0, llx = 0, lly = 0;
+  BBTopography_t *topo = NULL;
+  double RE = 0.0, R = 0;
+  PolarNavigator_t* navigator = NULL;
+  long bi = 0, ri = 0;
+  long nbins = 0, nrays = 0;
+  double* phi = NULL;
+  double* groundRange = NULL;
+  double height = 0.0;
+  double beamwidth = 0.0, elangle = 0.0;
+  double c, elLim, bb_tot;
+  double elBlock = 0.0;
 
   RAVE_ASSERT((self != NULL), "self == NULL");
+
   if (scan == NULL) {
     return NULL;
   }
+  navigator = PolarScan_getNavigator(scan);
+  if (navigator == NULL) {
+    RAVE_ERROR0("Scan does not have a polar navigator instance attached");
+    goto done;
+  }
 
-  topo = BeamBlockageMap_readTopography(self->mapper,
-                                        PolarScan_getLatitude(scan),
-                                        PolarScan_getLongitude(scan),
-                                        PolarScan_getMaxDistance(scan));
+  groundRange = BeamBlockageInternal_computeGroundRange(self, scan);
+  if (groundRange == NULL) {
+    goto done;
+  }
+
+  topo = BeamBlockageMap_getTopographyForScan(self->mapper, scan);
   if (topo == NULL) {
     goto done;
   }
 
-#ifdef KALLE
-  if (!BeamBlockageInternal_getLonLatFields(scan, &lon, &lat)) {
+  nbins = PolarScan_getNbins(scan);
+  nrays = PolarScan_getNrays(scan);
+
+  field = RAVE_OBJECT_NEW(&RaveField_TYPE);
+  if (field == NULL || !RaveField_createData(field, nbins, nrays, RaveDataType_DOUBLE)) {
     goto done;
   }
 
-  if (!BeamBlockageInternal_min(lon, &lon_min) ||
-      !BeamBlockageInternal_max(lon, &lon_max) ||
-      !BeamBlockageInternal_min(lat, &lat_min) ||
-      !BeamBlockageInternal_max(lon, &lat_max)) {
+  phi = RAVE_MALLOC(sizeof(double)*nbins);
+  if (phi == NULL) {
     goto done;
   }
-  lon_min = floor((lon_min - BBTopography_getUlxmap(topo)) / BBTopography_getXDim(topo)) - 1;
-  lon_max = ceil((lon_max - BBTopography_getUlxmap(topo)) / BBTopography_getXDim(topo)) + 1;
-  lat_max = floor((BBTopography_getUlymap(topo) - lat_max) / BBTopography_getYDim(topo)) - 1;
-  lat_min = ceil((BBTopography_getUlymap(topo) - lat_min) / BBTopography_getYDim(topo)) + 1;
-
-#endif
-  /*
-  lon_min = floor( (min_double(lon,ri*ai)-ulxmap) / xdim ) - 1;
-  lon_max = ceil( (max_double(lon,ri*ai)-ulxmap) / xdim ) + 1;
-  lat_max = floor( (ulymap-max_double(lat,ri*ai)) / ydim ) - 1;
-  lat_min = ceil( (ulymap-min_double(lat,ri*ai)) / ydim ) + 1;
-  */
+
+  RE = PolarNavigator_getEarthRadiusOrigin(navigator);
+  R = 1.0/((1.0/RE) + PolarNavigator_getDndh(navigator));
+  height = PolarNavigator_getAlt0(navigator);
+
+  beamwidth = PolarScan_getBeamwidth(scan);
+  elangle = PolarScan_getElangle(scan);
+
+  /* Width of Gaussian */
+  c = -((beamwidth/2.0)*(beamwidth/2.0))/log(0.5);
+
+  /* Elevation limits */
+  elLim = sqrt( -c*log(pow(10.0,(dBlim/10.0)) ) );
+
+  /* Find total blockage within -elLim to +elLim */
+  bb_tot = sqrt(M_PI*c) * erf(elLim/sqrt(c));
+
+  for (ri = 0; ri < nrays; ri++) {
+    for (bi = 0; bi < nbins; bi++) {
+      double v = 0.0;
+      BBTopography_getValue(topo, bi, ri, &v);
+      phi[bi] = asin((((v+R)*(v+R)) - (groundRange[bi]*groundRange[bi]) - ((R+height)*(R+height))) / (2*groundRange[bi]*(R+height)));
+      //fprintf(stderr, "ri=%d, bi=%d => phi(%d) = %f\n", ri, bi, bi, phi[bi]);
+      //phi[j] = asin((pow(v + R, 2) - pow(groundRange[bi],2) - pow(R + height,2) ) / (2*groundRange[bi]*(R+height)));
+    }
+    BeamBlockageInternal_cummax(phi, nbins);
+
+    for (bi = 0; bi < nbins; bi++) {
+      double bbval = 0.0;
+      elBlock = phi[bi];
+      if (elBlock < elangle - elLim) {
+        elBlock = -9999.0;
+      }
+      if (elBlock > elangle + elLim) {
+        elBlock = elangle + elLim;
+      }
+      bbval = -1.0/2.0 * sqrt(M_PI * c) * (erf((elangle - elBlock)/sqrt(c)) - erf(elLim/sqrt(c)))/bb_tot;
+      RaveField_setValue(field, bi, ri, bbval);
+      /*bb[i*ri+j] = -1./2. * sqrt(M_PI*c) * ( erf( (el-elBlock[i*ri+j]) / \
+                      sqrt(c) ) - erf( (el-(el-elLim))/sqrt(c) ) ) / bb_tot; */
+    }
+  }
+
   result = RAVE_OBJECT_COPY(field);
 done:
+  RAVE_OBJECT_RELEASE(navigator);
   RAVE_OBJECT_RELEASE(field);
-  RAVE_OBJECT_RELEASE(lon);
-  RAVE_OBJECT_RELEASE(lat);
+  RAVE_OBJECT_RELEASE(topo);
+  RAVE_FREE(phi);
+  RAVE_FREE(groundRange);
   return result;
 }
 
index 361891f..37b89a9 100644 (file)
@@ -59,8 +59,9 @@ const char* BeamBlockage_getTopo30Directory(BeamBlockage_t* self);
  * Gets the blockage for the provided scan.
  * @param[in] self - self
  * @param[in] scan - the scan to check blockage
+ * @param[in] dBlim -
  * @return the beam blockage field
  */
-RaveField_t* BeamBlockage_getBlockage(BeamBlockage_t* self, PolarScan_t* scan);
+RaveField_t* BeamBlockage_getBlockage(BeamBlockage_t* self, PolarScan_t* scan, double dBlim);
 
 #endif /* BEAMBLOCKAGE_H */
index 01e9d7f..6ba1910 100644 (file)
@@ -152,13 +152,13 @@ static BBTopography_t* BeamBlockageMapInternal_readHeader(BeamBlockageMap_t* sel
       } else if (strcmp(token, "NBITS") == 0) {
         nbits = (int)f;
       } else if (strcmp(token, "ULXMAP") == 0) {
-        BBTopography_setUlxmap(field, f);
+        BBTopography_setUlxmap(field, f * M_PI/180.0);
       } else if (strcmp(token, "ULYMAP") == 0) {
-        BBTopography_setUlymap(field, f);
+        BBTopography_setUlymap(field, f * M_PI/180.0);
       } else if (strcmp(token, "XDIM") == 0) {
-        BBTopography_setXDim(field, f);
+        BBTopography_setXDim(field, f * M_PI/180.0);
       } else if (strcmp(token, "YDIM") == 0) {
-        BBTopography_setYDim(field, f);
+        BBTopography_setYDim(field, f * M_PI/180.0);
       }
     }
   }
@@ -184,6 +184,15 @@ done:
   return result;
 }
 
+/**
+ * Fills the data in the topography field. The topography field should first have been created
+ * with a call to \ref BeamBlockageMapInternal_readHeader in order to get correct dimensions
+ * etc.
+ * @param[in] self - self
+ * @param[in] filename - the filename (exluding suffix and directory name)
+ * @param[in] field - the gtopo30 topography skeleton.
+ * @returns 1 on success otherwise 0
+ */
 static int BeamBlockageMapInternal_fillData(BeamBlockageMap_t* self, const char* filename, BBTopography_t* field)
 {
   int result = 0;
@@ -245,6 +254,13 @@ done:
   return result;
 }
 
+/**
+ * Reads the topography with the specified filename. The directory path should not be used, instead
+ * the topo30dir variable will be used.
+ * @param[in] self - self
+ * @param[in] filename - the gtopo30 filename (excluding suffix)
+ * @returns the topography on success otherwise NULL
+ */
 static BBTopography_t* BeamBlockageMapInternal_readTopography(BeamBlockageMap_t* self, const char* filename)
 {
   BBTopography_t *result = NULL, *field = NULL;
@@ -265,149 +281,33 @@ done:
   RAVE_OBJECT_RELEASE(field);
   return result;
 }
-/**
- * Locates the min value in the rave data2d field.
- * @param[in] field - the data field
- * @param[out] omin - the minimum value
- * @return 1 on success otherwise 0
- */
-static int BeamBlockageMapInternal_min(RaveData2D_t* field, double* omin)
-{
-  double minv = 0.0, v = 0.0;
-  long xsize = 0, ysize = 0;
-  long x = 0, y = 0;
-  int result = 0;
-
-  RAVE_ASSERT((field != NULL), "field == NULL");
-  RAVE_ASSERT((omin != NULL), "omin == NULL");
-
-  xsize = RaveData2D_getXsize(field);
-  ysize = RaveData2D_getYsize(field);
-  if (x > 0 && y > 0) {
-    result = RaveData2D_getValue(field, 0, 0, &minv);
-    for (x = 0; result && x < xsize; x++) {
-      for (y = 0; result && y < ysize; y++) {
-        result = RaveData2D_getValue(field, x, y, &v);
-        if (v <= minv) {
-          minv = v;
-        }
-      }
-    }
-  }
-  *omin = minv;
-  return result;
-}
-
-/**
- * Locates the max value in the rave data2d field.
- * @param[in] field - the data field
- * @param[out] omax - the minimum value
- * @return 1 on success otherwise 0
- */
-static int BeamBlockageMapInternal_max(RaveData2D_t* field, double* omax)
-{
-  double maxv = 0.0, v = 0.0;
-  long xsize = 0, ysize = 0;
-  long x = 0, y = 0;
-  int result = 0;
-
-  RAVE_ASSERT((field != NULL), "field == NULL");
-  RAVE_ASSERT((omax != NULL), "omax == NULL");
-
-  xsize = RaveData2D_getXsize(field);
-  ysize = RaveData2D_getYsize(field);
-  if (x > 0 && y > 0) {
-    result = RaveData2D_getValue(field, 0, 0, &maxv);
-    for (x = 0; result && x < xsize; x++) {
-      for (y = 0; result && y < ysize; y++) {
-        result = RaveData2D_getValue(field, x, y, &v);
-        if (v >= maxv) {
-          maxv = v;
-        }
-      }
-    }
-  }
-  *omax = maxv;
-  return result;
-}
 
 /**
- * Extracts a subset of a rave data 2d field
- * @param[in] field - the original field
- * @param[in] ulx - the upper left x coordinate
- * @param[in] uly - the upper left y coordinate
- * @param[in] llx - the lower left x coordinate
- * @param[in] lly - the lower left y coordinate
- * @returns the subset on success otherwise NULL
+ * Creates a topography that is mapped against a specific scan.
+ * @param[in] self - self
+ * @param[in] topo - the overall topography that hopefully covers the scan
+ * @param[in] scan - the scan that should get the topography mapped
+ * @return the mapped topography on success otherwise NULL
  */
-static RaveData2D_t* BeamBlockageMapInternal_subset(RaveData2D_t* field, long ulx, long uly, long llx, long lly)
-{
-  RaveData2D_t *ofield = NULL, *result = NULL;
-  long xsize = 0, ysize = 0, x = 0, y = 0;
-  RaveDataType dtype = RaveDataType_UNDEFINED;
-
-  RAVE_ASSERT((field != NULL), "field == NULL");
-
-  xsize = llx - ulx;
-  ysize = lly - uly;
-  dtype = RaveData2D_getType(field);
-
-  if (ulx < 0 || uly < 0 ||
-      llx >= RaveData2D_getXsize(field) ||
-      lly >= RaveData2D_getYsize(field) ||
-      xsize <= 0 || ysize <= 0) {
-    RAVE_WARNING0("Bad subset coorner points.");
-    goto done;
-  }
-
-  ofield = RAVE_OBJECT_NEW(&RaveData2D_TYPE);
-  if (ofield == NULL || !RaveData2D_createData(ofield, xsize, ysize, dtype)) {
-    goto done;
-  }
-
-  for (x = ulx; x < llx; x++) {
-    for (y = uly; y < lly; y++) {
-      double v = 0.0;
-      RaveData2D_getValue(field, x, y, &v);
-      RaveData2D_setValue(ofield, x-ulx, y-uly, v);
-    }
-  }
-
-  result = RAVE_OBJECT_COPY(ofield);
-done:
-  RAVE_OBJECT_RELEASE(ofield);
-  return result;
-}
-
-
-static int BeamBlockageMapInternal_getLonLatFields(PolarScan_t* scan, RaveData2D_t** lon, RaveData2D_t** lat)
+static BBTopography_t* BeamBlockageMapInternal_createMappedTopography(BeamBlockageMap_t* self, BBTopography_t* topo, PolarScan_t* scan)
 {
-  int result = 0;
+  BBTopography_t *field = NULL, *result = NULL;
   long nrays = 0, nbins = 0;
   long ri = 0, bi = 0;
-  double elangle = 0.0;
-  RaveData2D_t *olon = NULL, *olat = NULL;
 
-  RAVE_ASSERT((lon != NULL), "lon == NULL");
-  RAVE_ASSERT((lat != NULL), "lat == NULL");
-
-  if (scan == NULL) {
-    goto done;
-  }
+  RAVE_ASSERT((self != NULL), "self == NULL");
+  RAVE_ASSERT((topo != NULL), "topo == NULL");
+  RAVE_ASSERT((scan != NULL), "scan == NULL");
 
   nrays = PolarScan_getNrays(scan);
   nbins = PolarScan_getNbins(scan);
-  elangle = PolarScan_getElangle(scan);
-
-  olon = RAVE_OBJECT_NEW(&RaveData2D_TYPE);
-  olat = RAVE_OBJECT_NEW(&RaveData2D_TYPE);
 
-  if (olon == NULL || olat == NULL) {
+  field = RAVE_OBJECT_NEW(&BBTopography_TYPE);
+  if (field == NULL) {
     goto done;
   }
-
-  if (!RaveData2D_createData(olon, nbins, nrays, RaveDataType_DOUBLE) ||
-      !RaveData2D_createData(olat, nbins, nrays, RaveDataType_DOUBLE)) {
+  if (!BBTopography_createData(field, nbins, nrays, BBTopography_getDataType(topo))) {
+    RAVE_ERROR0("Failed to create data field");
     goto done;
   }
 
@@ -415,23 +315,23 @@ static int BeamBlockageMapInternal_getLonLatFields(PolarScan_t* scan, RaveData2D
     for (bi = 0; bi < nbins; bi++) {
       double lonval = 0.0, latval = 0.0;
       if (PolarScan_getLonLatFromIndex(scan, bi, ri, &lonval, &latval)) {
-        RaveData2D_setValue(olon, bi, ri, lonval);
-        RaveData2D_setValue(olat, bi, ri, latval);
+        double v = 0.0;
+        BBTopography_getValueAtLonLat(topo, lonval, latval, &v);
+        BBTopography_setValue(field, bi, ri, v);
       }
     }
   }
 
-  *lat = RAVE_OBJECT_COPY(olat);
-  *lon = RAVE_OBJECT_COPY(olon);
-
-  result = 1;
+  result = RAVE_OBJECT_COPY(field);
 done:
-  RAVE_OBJECT_RELEASE(olon);
-  RAVE_OBJECT_RELEASE(olat);
-
+  RAVE_OBJECT_RELEASE(field);
   return result;
 }
 
+/*@} End of Private functions */
+
+/*@{ Interface functions */
+
 /**
  * Find out which maps are needed to cover given area
  * @param[in] lat - latitude of radar in degrees
@@ -492,11 +392,7 @@ done:
 BBTopography_t* BeamBlockageMap_getTopographyForScan(BeamBlockageMap_t* self, PolarScan_t* scan)
 {
   BBTopography_t *field = NULL, *result = NULL;
-  RaveData2D_t *lon = NULL, *lat = NULL;
   BBTopography_t* topo = NULL;
-  double dlon_min = 0.0, dlon_max = 0.0, dlat_min = 0.0, dlat_max = 0.0;
-  double ddlon = 0.0, ddlat = 0.0;
-  long ulx = 0, uly = 0, llx = 0, lly = 0;
 
   RAVE_ASSERT((self != NULL), "self == NULL");
   if (scan == NULL) {
@@ -511,36 +407,15 @@ BBTopography_t* BeamBlockageMap_getTopographyForScan(BeamBlockageMap_t* self, Po
     goto done;
   }
 
-  if (!BeamBlockageMapInternal_getLonLatFields(scan, &lon, &lat)) {
-    goto done;
-  }
+  field = BeamBlockageMapInternal_createMappedTopography(self, topo, scan);
 
-  if (!BeamBlockageMapInternal_min(lon, &dlon_min) ||
-      !BeamBlockageMapInternal_max(lon, &dlon_max) ||
-      !BeamBlockageMapInternal_min(lat, &dlat_min) ||
-      !BeamBlockageMapInternal_max(lon, &dlat_max)) {
-    goto done;
-  }
-
-  ulx = floor((dlon_min - BBTopography_getUlxmap(topo)) / BBTopography_getXDim(topo)) - 1; /* Min lon */
-  uly = floor((BBTopography_getUlymap(topo) - dlat_max) / BBTopography_getYDim(topo)) - 1; /* Max lat */
-
-  llx = ceil((dlon_max - BBTopography_getUlxmap(topo)) / BBTopography_getYDim(topo)) + 1;  /* Max Lon */
-  //lly = ceil(())
-
-  /*
-  lon_min = floor((dlon_min - BBTopography_getUlxmap(topo)) / BBTopography_getXDim(topo)) - 1;
-  lon_max = ceil((dlon_max - BBTopography_getUlxmap(topo)) / BBTopography_getXDim(topo)) + 1;
-  lat_max = floor((BBTopography_getUlymap(topo) - dlat_max) / BBTopography_getYDim(topo)) - 1;
-  lat_min = ceil((BBTopography_getUlymap(topo) - dlat_min) / BBTopography_getYDim(topo)) + 1;
-  */
   result = RAVE_OBJECT_COPY(field);
 done:
   RAVE_OBJECT_RELEASE(field);
+  RAVE_OBJECT_RELEASE(topo);
   return result;
 }
 
-/*@{ Interface functions */
 int BeamBlockageMap_setTopo30Directory(BeamBlockageMap_t* self, const char* topodirectory)
 {
   char* tmp = NULL;
index 1c5d092..d6e280a 100644 (file)
@@ -30,6 +30,8 @@ along with beamb.  If not, see <http://www.gnu.org/licenses/>.
 #define PYBEAMBLOCKAGE_MODULE   /**< to get correct part in pybeamblockage */
 #include "pybeamblockage.h"
 
+#include "pypolarscan.h"
+#include "pyravefield.h"
 #include "pyrave_debug.h"
 #include "rave_alloc.h"
 
@@ -144,11 +146,41 @@ static PyObject* _pybeamblockage_new(PyObject* self, PyObject* args)
 }
 
 /**
+ * Returns the blockage for the provided scan given gaussian limit.
+ * @param[in] self - self
+ * @param[in] args - the arguments (PyPolarScan, double (Limit of Gaussian approximation of main lobe))
+ * @return the PyRaveField on success otherwise NULL
+ */
+static PyObject* _pybeamblockage_getBlockage(PyBeamBlockage* self, PyObject* args)
+{
+  PyObject* pyin = NULL;
+  double dBlim = 0;
+  RaveField_t* field = NULL;
+  PyObject* result = NULL;
+
+  if (!PyArg_ParseTuple(args, "Od", &pyin, &dBlim)) {
+    return NULL;
+  }
+
+  if (!PyPolarScan_Check(pyin)) {
+    raiseException_returnNULL(PyExc_ValueError, "First argument should be a Polar Scan");
+  }
+
+  field = BeamBlockage_getBlockage(self->beamb, ((PyPolarScan*)pyin)->scan, dBlim);
+  if (field != NULL) {
+    result = (PyObject*)PyRaveField_New(field);
+  }
+  RAVE_OBJECT_RELEASE(field);
+  return result;
+}
+
+/**
  * All methods a ropo generator can have
  */
 static struct PyMethodDef _pybeamblockage_methods[] =
 {
   {"topo30dir", NULL},
+  {"getBlockage", (PyCFunction)_pybeamblockage_getBlockage, 1},
   {NULL, NULL} /* sentinel */
 };
 
@@ -270,6 +302,8 @@ init_beamblockage(void)
   if (ErrorObject == NULL || PyDict_SetItemString(dictionary, "error", ErrorObject) != 0) {
     Py_FatalError("Can't define _beamblockage.error");
   }
+  import_pyravefield();
+  import_pypolarscan();
   PYRAVE_DEBUG_INITIALIZE;
 }
 /*@} End of Module setup */
index 2420ce2..5955d3e 100644 (file)
@@ -30,6 +30,7 @@ along with beamb.  If not, see <http://www.gnu.org/licenses/>.
 #define PYBEAMBLOCKAGEMAP_MODULE   /**< to get correct part in pybeamblockagemap */
 #include "pybeamblockagemap.h"
 
+#include "pypolarscan.h"
 #include "pyrave_debug.h"
 #include "rave_alloc.h"
 #include "pybbtopography.h"
@@ -170,6 +171,25 @@ static PyObject* _pybeamblockagemap_readTopography(PyBeamBlockageMap* self, PyOb
   return result;
 }
 
+static PyObject* _pybeamblockagemap_getTopographyForScan(PyBeamBlockageMap* self, PyObject* args)
+{
+  PyObject* pyin = NULL;
+  BBTopography_t* field = NULL;
+  PyObject* result = NULL;
+  if (!PyArg_ParseTuple(args, "O", &pyin)) {
+    return NULL;
+  }
+  if (!PyPolarScan_Check(pyin)) {
+    raiseException_returnNULL(PyExc_ValueError, "In object must be a polar scan");
+  }
+  field = BeamBlockageMap_getTopographyForScan(self->map, ((PyPolarScan*)pyin)->scan);
+  if (field != NULL) {
+    result = (PyObject*)PyBBTopography_New(field);
+  }
+  RAVE_OBJECT_RELEASE(field);
+  return result;
+}
+
 /**
  * All methods a beam blockage map can have
  */
@@ -177,6 +197,7 @@ static struct PyMethodDef _pybeamblockagemap_methods[] =
 {
   {"topo30dir", NULL},
   {"readTopography", (PyCFunction)_pybeamblockagemap_readTopography, 1},
+  {"getTopographyForScan", (PyCFunction)_pybeamblockagemap_getTopographyForScan, 1},
   {NULL, NULL} /* sentinel */
 };
 
@@ -299,6 +320,7 @@ init_beamblockagemap(void)
     Py_FatalError("Can't define _beamblockagemap.error");
   }
   import_bbtopography();
+  import_pypolarscan();
   PYRAVE_DEBUG_INITIALIZE;
 }
 /*@} End of Module setup */
index 34ff32c..99bbd6c 100644 (file)
@@ -33,6 +33,8 @@ import _rave
 import math
 
 class PyBeamBlockageMapTest(unittest.TestCase):
+  SCAN_FILENAME = "fixtures/scan_sevil_20100702T113200Z.h5"
+  
   def setUp(self):
     pass    
 
@@ -65,6 +67,13 @@ class PyBeamBlockageMapTest(unittest.TestCase):
     self.assertEquals(6000, result.nrows)
     self.assertEquals(9600, result.ncols)
   
+  def testGetTopographyForScan(self):
+    a = _beamblockagemap.new()
+    a.topo30dir="../../data/gtopo30"
+    scan = _raveio.open(self.SCAN_FILENAME).object
+    topo = a.getTopographyForScan(scan)
+    self.assertEqual(174, topo.getData()[0][0])
+    
 if __name__ == "__main__":
   #import sys;sys.argv = ['', 'Test.testName']
   unittest.main()
\ No newline at end of file
index 35eef21..52c0ed7 100644 (file)
@@ -31,6 +31,7 @@ import os, string
 import _rave
 
 class PyBeamBlockageTest(unittest.TestCase):
+  SCAN_FILENAME = "fixtures/scan_sevil_20100702T113200Z.h5"
   def setUp(self):
     pass    
 
@@ -51,6 +52,13 @@ class PyBeamBlockageTest(unittest.TestCase):
     a = _beamblockage.new()
     a.topo30dir="../../data/gtopo30"
     
+  def test_getBlockage(self):
+    a = _beamblockage.new()
+    a.topo30dir="../../data/gtopo30"
+    scan = _raveio.open(self.SCAN_FILENAME).object
+    result = a.getBlockage(scan, -20.0);
+    print `result.getData()`
+    
   
 if __name__ == "__main__":
   #import sys;sys.argv = ['', 'Test.testName']
diff --git a/test/pytest/fixtures/scan_sevil_20100702T113200Z.h5 b/test/pytest/fixtures/scan_sevil_20100702T113200Z.h5
new file mode 100644 (file)
index 0000000..71baffb
Binary files /dev/null and b/test/pytest/fixtures/scan_sevil_20100702T113200Z.h5 differ