A code block in the cliplib was surrounded with "try" and "finally" to make...

A code block in the cliplib was surrounded with "try" and "finally" to make sure the raster is always reset whatever happens.

toolbox/cliplib.py

@@ -132,32 +132,34 @@ def get_masked_clip_result(datagrid, nvlp, mask, initvalue=0):
     # Extract data from the datagrid and return it as a 2-dimensional array
     shape = (nvlp.nrows, nvlp.ncols)
     datatype = get_datatype(datagrid.dataformat)
-    if initvalue == 0: result = np.zeros(shape, dtype=datatype)
+    result = np.zeros(shape, dtype=datatype)
+    if initvalue != 0: result = initvalue
     
-    # Now calculate which part of the datagrid we need
-    # Assume that the datagrid coincides exactly with the "country pixels" 
-    mink = int(round((nvlp.xll - datagrid.xll) / nvlp.dx))
-    maxk = mink + nvlp.ncols - 1
-    maxi = datagrid.nrows - int(round((nvlp.yll - datagrid.yll) / nvlp.dy)) - 1
-    mini = maxi - nvlp.nrows + 1
-    
-    # Loop over the rows - assume that the file position is zero
-    for i in range(datagrid.nrows):
-        if i < mini: 
-            datagrid.next(False)
-            continue
-        if i > maxi: break
-        line = datagrid.next()
-        arrline = line[mink:maxk+1]
-        if isinstance(arrline, array):
-            arrline = np.array(arrline) 
-        arrline.shape = (nvlp.ncols)
-        result[i - mini, :] = arrline
+    try:
+        # Now calculate which part of the datagrid we need
+        # Assume that the datagrid coincides exactly with the "country pixels" 
+        mink = int(round((nvlp.xll - datagrid.xll) / nvlp.dx))
+        maxk = mink + nvlp.ncols - 1
+        maxi = datagrid.nrows - int(round((nvlp.yll - datagrid.yll) / nvlp.dy)) - 1
+        mini = maxi - nvlp.nrows + 1
         
-    # Wind up
-    datagrid.reset()
-    result = ma.masked_array(result, mask)
-    return result
+        # Loop over the rows - assume that the file position is zero
+        for i in range(datagrid.nrows):
+            if i < mini: 
+                datagrid.next(False)
+                continue
+            if i > maxi: break
+            line = datagrid.next()
+            arrline = line[mink:maxk+1]
+            if isinstance(arrline, array):
+                arrline = np.array(arrline) 
+            arrline.shape = (nvlp.ncols)
+            result[i - mini, :] = arrline
+    finally:
+        # Wind up
+        datagrid.reset()
+        result = ma.masked_array(result, mask)
+        return result
         
 def most_left(line, value):
     result = -1
@@ -238,7 +240,7 @@ def join_envelopes(nvlps):
     return result
            
 # Method create_raster that is to return an envelope and an array with data
-# Input is a minimum envelope and a masked array with values
+# Input is a minimal envelope and a masked array with values
 def create_raster(countrydata, envelope=None, nodatavalue=None, debug="n"):
     # Assume that the countrydata are tuples with envelope and masked arrays
     if envelope == None:

toolbox/punchlib.py

+import numpy as np
+import lmgeo.toolbox.cliplib as cp
+import lmgeo.formats.gridenvelope2d as ge
+from pyproj import Proj, transform
+
+# Constants
+eps = 0.00000001
+
+# Determine if a point is inside a given polygon or not
+# Polygon is a list of (x,y) pairs. See http://www.ariel.com.au/a/python-point-int-poly.html
+def point_inside_polygon(x,y,poly):
+    n = len(poly)
+    inside =False
+
+    p1x,p1y = poly[0]
+    for i in range(n+1):
+        p2x,p2y = poly[i % n]
+        if y > min(p1y,p2y):
+            if y <= max(p1y,p2y):
+                if x <= max(p1x,p2x):
+                    if p1y != p2y:
+                        xinters = (y-p1y)*(p2x-p1x)/(p2y-p1y)+p1x
+                    if p1x == p2x or x <= xinters:
+                        inside = not inside
+        p1x,p1y = p2x,p2y
+    return inside
+
+def point_inside_tetragon(x,y,poly):
+    # Initialise
+    result = False
+    
+    # Determine twice two lines and check that the point lies in between them
+    # TODO: the tetragon could be self-intersecting - check for this
+    for cornerpairs in [[(0,1), (2,3)], [(1,2), (3,0)]]:
+        params = []
+        for (i, j) in cornerpairs:
+            params.append(get_linear_coeffs(poly[i], poly[j]))
+        result = point_between_lines(x, y, params[0], params[1])
+        if not result: break
+    return result      
+
+def point_between_lines(x, y, params1, params2):
+    if params1[1] == float("inf") or params2[1] == float("inf"):
+        params1 = swap_axes(params1)
+        params2 = swap_axes(params2)
+        x, y = y, x
+    y1 = params1[0] + params1[1] * x 
+    y2 = params2[0] + params2[1] * x
+    result = ((y1 <= y) and (y2 >= y)) or ((y1 >= y) and (y2 <= y)) 
+    return result
+
+def swap_axes(params):
+    if params[1] == float("inf"): params[1] = 0.0
+    else: 
+        if abs(params[1]) > eps: params[1] = 1 / float(params[1])
+        else: params[1] = 1 / eps
+        params[0] = -1 * params[0] * params[1]
+    return params
+    
+def get_linear_coeffs(pt1, pt2):
+    dx = float(pt2[0] - pt1[0])
+    if abs(dx) > eps: 
+        rc = (pt2[1] - pt1[1]) / dx 
+        ic = pt2[1] - rc * pt2[0]
+    else: 
+        rc = float("inf")
+        ic = pt2[0]
+    return [ic, rc] 
+
+def snap_feature_bounds(bbox, rastergrid):
+    # Initialise
+    result = 4 * [0.0]
+
+    # Now calculate the offsets - bbox has to be in the same coordinate system!
+    factor = divmod(bbox[0] - rastergrid.xll, rastergrid.dx)[0]
+    result[0] = rastergrid.xll + factor * rastergrid.dx
+    factor = divmod(bbox[1] - rastergrid.yll, rastergrid.dy)[0]
+    result[1] = rastergrid.yll + factor * rastergrid.dy
+    factor = 1.0 + divmod(bbox[2] - rastergrid.xll, rastergrid.dx)[0]
+    result[2] = rastergrid.xll + factor * rastergrid.dx
+    factor = 1.0 + divmod(bbox[3] - rastergrid.yll, rastergrid.dy)[0]
+    result[3] = rastergrid.yll + factor * rastergrid.dy
+    return result
+
+"""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]
+    dx, dy = 2 * [cellsize]
+    
+    # Initialise further
+    ncols = int((snapped_bbox[2] - xll) / dx)
+    nrows = int((snapped_bbox[3] - yll) / dy)
+    result = np.zeros((nrows, ncols), dtype=np.int)
+    
+    try:      
+        # Determine the nature of the polygon
+        tetragon = False
+        if len(points) == 4:
+            tetragon = True
+        elif len(points) == 5:
+            pt1, ptn = points[0], points[-1]
+            if (abs(pt1[0] - ptn[0]) < eps) and (abs(pt1[1] - ptn[1]) < eps):
+                tetragon = True            
+        
+        # Now pass over the grid within the boundary box - punch where centroids are located
+        for i in range(nrows):
+            for k in range(ncols):
+                centroid = (xll + (k+0.5)*dx, yll + (nrows-i-0.5)*dy)
+                if tetragon:
+                    result[i, k] = int(not point_inside_tetragon(centroid[0], centroid[1], points))
+                else:
+                    result[i, k] = int(not point_inside_polygon(centroid[0], centroid[1], points))
+    except Exception as e:    
+        print str(e)        
+    finally:
+        return result
+
+"""This is to find the index of a field in the fields of a shapefile or similar structure"""
+def get_field_index(fields, srchfld):
+    # Check whether the searched field is in the fields
+    k = 0
+    for field in fields:
+        if field[0] == srchfld:
+            found = True
+            break
+        k += 1
+    if not found:
+        raise ValueError("Indicated field %s is not part of the shapefile" % srchfld)
+    k = k - 1 # DeletionFlag is always the first field, but is not represented in the records
+    return k
+
+"""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 
+   values are not interesting, then f could be: lambda x: divmod(x, 10)[0] """
+def get_zonal_stats_record(myshape, rastergrid, stat_types, inproj4='', outproj4='', f=None, idname=''):
+    # Initialise and check input
+    if "TABLE" in stat_types:
+        result = []
+        if len(stat_types) > 1:
+            msg = "A table of values and their counts cannot be generated together with other statistics."
+            raise Warning(msg)
+            stat_types = ["TABLE"]
+    else:
+        result = len(stat_types) * [0.0]  
+    
+    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
+            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])        
+        else:
+            bbox = myshape.bbox
+            points = myshape.points            
+        
+        # Now get an area of whole cells that completely cover this area
+        snapped_bbox = snap_feature_bounds(bbox, rastergrid)
+        
+        # Mask the rastergrid
+        mask = get_mask(snapped_bbox, points, rastergrid.cellsize, inproj4, outproj4)
+        (nrows, ncols) = mask.shape
+        (xll, yll) = snapped_bbox[0:2]
+        nvlp = ge.GridEnvelope2D(ncols, nrows, xll, yll, rastergrid.cellsize, rastergrid.cellsize)  
+        masked_data = cp.get_masked_clip_result(rastergrid, nvlp, mask)
+        assert masked_data.size > 0, "Array with masked data has zero size!"
+        
+        # Now obtain the statistics
+        stat = ""
+        for stat in stat_types:
+            value = 0.0
+            if stat == "MEAN":
+                value = np.mean(masked_data)
+            elif stat == "MAX":
+                value = np.max(masked_data)
+            elif stat == "MIN":
+                value = np.min(masked_data)
+            elif stat == "MEDIAN":
+                value = np.median(masked_data)
+            elif stat == "SUM":
+                value = np.sum(masked_data)
+            elif stat == "STD":
+                value = np.std(masked_data)
+            elif stat == "TABLE":
+                break
+            else:
+                raise ValueError("Indicated statistic not supported: " + stat)
+            
+        if stat != "TABLE":
+            result.append(value)
+        else:
+            # Retrieve the values and tally them (Dutch: turven)
+            # The array masked_data is a numpy masked array
+            if idname == '': idname = "value"
+            if f != None: masked_data = f(masked_data) 
+            for i in range(nrows):
+                for k in range(ncols):
+                    if mask[i,k] == 1: continue
+                    value = masked_data[i,k]
+                    valdict = filter(lambda x: x[idname] == value, result)
+                    if len(valdict) != 0: valdict[0]["count"] += 1
+                    else: result.append({idname: value, "count":1}) 
+            result = sorted(result, key=lambda k: k["count"]) 
+            
+    except Exception as e:
+        print str(e)
+    finally:
+        return result
+    
+def get_zonal_stats_as_table(shapeRecords, srchfld, rastergrid, stat_types):
+    result = []
+    rec = {srchfld:0}
+    k = get_field_index(shapeRecords[0], srchfld)
+                    
+    for shpRec in shapeRecords:
+        fldvalue = shpRec.record[k]
+        statrec = get_zonal_stats_record(shpRec.shape, rastergrid, stat_types)
+        rec[srchfld] = fldvalue
+        for i in range(len(stat_types)):
+            stat = stat_types[i]
+            rec[stat] = statrec[i]
+        result.append(rec)
+    return result
+
+
+        
+    
\ No newline at end of file