method get_transformed_boundary_box defined separately

toolbox/punchlib.py

-import numpy as np 
+import numpy as np
 import lmgeo.toolbox.cliplib as cp
 import lmgeo.formats.gridenvelope2d as ge
 from pyproj import Proj, transform
@@ -82,7 +82,7 @@ def snap_feature_bounds(bbox, rastergrid):
     result[3] = rastergrid.yll + factor * rastergrid.dy
     return result
 
-"""Avoid calling the following method for too large boundary boxes to avoid slow performance"""
+"""Avoid calling the following method for too large boundary boxes, to avoid slow performance"""
 def get_mask(snapped_bbox, points, cellsize, inproj4='', outproj4=''):
     # Assume that the cellsize can be derived from the rastergrid
     xll, yll = snapped_bbox[0], snapped_bbox[1]
@@ -130,6 +130,41 @@ def get_field_index(fields, srchfld):
     k = k - 1 # DeletionFlag is always the first field, but is not represented in the records
     return k
 
+def get_transformed_boundary_box(bbox, inproj4='', outproj4='', rotate=True):
+    # Initialise
+    result = []
+    
+    # Check input
+    if (inproj4 == '') or (outproj4 == ''):
+        raise Exception("Proj.4 string cannot be empty!")
+    
+    try:
+        inProj = Proj(inproj4)
+        outProj = Proj(outproj4)
+        pt1 = (bbox[0], bbox[1]) # LL
+        pt2 = (bbox[0], bbox[3]) # UL
+        pt3 = (bbox[2], bbox[1]) # LR
+        pt4 = (bbox[2], bbox[3]) # UR
+        minx, miny = transform(inProj, outProj, pt1[0], pt1[1])
+        maxx, maxy = minx, miny
+        if not rotate: result.append((minx,miny))
+        for pt in [pt2, pt3, pt4]:
+            x, y = transform(inProj, outProj, pt[0], pt[1])
+            if rotate:
+                minx = min(minx, x)
+                maxx = max(maxx, x)
+                miny = min(miny, y)
+                maxy = max(maxy, y) 
+            else:
+                result.append((x,y))
+        if rotate:         
+            result = [minx, miny, maxx, maxy] 
+                
+    except Exception as e:
+        print "Error in method get_transformed_boundary_box of module punchlib: " + str(e)
+    finally:
+        return result
+
 """Extra functionality is that the shape and the raster can have different coordinate systems
    and the obtained values can be transformed using the supplied function f before analysis,
    i.e. in case a table is requested showing a frequency analysis. E.g. if the last digit of the 
@@ -147,24 +182,12 @@ def get_zonal_stats_record(myshape, rastergrid, stat_types, inproj4='', outproj4
     
     try:
         if (inproj4 != '') and (outproj4 != ''):
-            # Apparently, the coordinate system is not the same for the raster and the grid
-            # Convert all the points in order to be sure we get a proper boundary box and mask
+            # Apparently, the coordinate system is not the same for the raster and the shape
+            # Convert all the points in order to be sure we get a proper boundary box and mask    
+            bbox = get_transformed_boundary_box(myshape.bbox)
+            points = []
             inProj = Proj(inproj4)
             outProj = Proj(outproj4)
-            pt1 = (myshape.bbox[0], myshape.bbox[1])
-            pt2 = (myshape.bbox[0], myshape.bbox[3])
-            pt3 = (myshape.bbox[2], myshape.bbox[1])
-            pt4 = (myshape.bbox[2], myshape.bbox[3])
-            minx, miny = transform(inProj, outProj, pt1[0], pt1[1])
-            maxx, maxy = minx, miny
-            for pt in [pt2, pt3, pt4]:
-                x, y = transform(inProj, outProj, pt[0], pt[1])
-                minx = min(minx, x)
-                maxx = max(maxx, x)
-                miny = min(miny, y)
-                maxy = max(maxy, y)       
-            bbox = [minx, miny, maxx, maxy] 
-            points = []
             for pt in myshape.points:
                 x, y = transform(inProj, outProj, pt[0], pt[1])
                 points.append([x, y])        
@@ -224,7 +247,8 @@ def get_zonal_stats_record(myshape, rastergrid, stat_types, inproj4='', outproj4
         print str(e)
     finally:
         return result
-    
+
+"""For each shape we'll tally how many grid cells with values a, b etc. occur within its boundary box"""
 def get_zonal_stats_as_table(shapeRecords, srchfld, rastergrid, stat_types):
     result = []
     rec = {srchfld:0}
@@ -240,6 +264,9 @@ def get_zonal_stats_as_table(shapeRecords, srchfld, rastergrid, stat_types):
         result.append(rec)
     return result
 
-
+def get_centroid(rastergrid, i, k):
+    result = [rastergrid.xll + (k+0.5)*rastergrid.dx]
+    result.extend([rastergrid.yll + (rastergrid.nrows-i-0.5)*rastergrid.dy])
+    return result
         
     
\ No newline at end of file