Merge pull request #167 from Hydrata/master

allow trim of internal_holes (houses) when producing geotiff outputs using Make_Geotif

@davekennewell  it seems that appveyor is playing up or at least conda-forge updating isn't working. This seems independent of anuga and your change. Your changes are passing travis. So I will merge your PR

anuga/geospatial_data/tests/test_geospatial_data.py

@@ -33,10 +33,12 @@ def test_0(self):
 
         # Check __repr__
         # FIXME (Ole): Is this really machine independent?
-        rep = `G`
+        rep = G.__repr__()
         ref = '[[ 1.   2.1]\n [ 3.   5.3]]'
-        msg = 'Representation %s is not equal to %s' % (rep, ref)
-        assert rep == ref, msg
+        # it seems the leading spaces are removed on some conda builds, so lets allow that too:
+        ref_conda = '[[1.  2.1]\n [3.  5.3]]'
+        msg = 'Representation %s is not equal to %s or %s' % (rep, ref, ref_conda)
+        assert ((rep == ref) or (rep == ref_conda)), msg
 
         #Check getter
         assert num.allclose(G.get_data_points(), [[1.0, 2.1], [3.0, 5.3]])

anuga/utilities/plot_utils.py

@@ -957,6 +957,7 @@ def Make_Geotif(swwFile=None,
              min_allowed_height=1.0e-05,
              output_dir='TIFS',
              bounding_polygon=None,
+             internal_holes=None,
              verbose=False,
              k_nearest_neighbours=3,
              creation_options=[]):
@@ -983,6 +984,7 @@ def Make_Geotif(swwFile=None,
                 min_allowed_height -- Minimum allowed height from ANUGA
                 output_dir -- Write outputs to this directory
                 bounding_polygon -- polygon (e.g. from read_polygon) If present, only set values of raster cells inside the bounding_polygon
+                internal_holes -- a list of polygons. If present, do not set values of raster cells inside these polygons.
                 k_nearest_neighbours -- how many neighbours to use in interpolation. If k>1, inverse-distance-weighted interpolation is used
                 creation_options -- list of tif creation options for gdal, e.g. ['COMPRESS=DEFLATE']
     """
@@ -1104,11 +1106,19 @@ def myInterpFun(quantity):
                 num += quantity[nn_inds[:,i]]*nn_wts[:,i]
             return (num/denom)
 
-    if(bounding_polygon is not None):
+    if bounding_polygon is not None:
         # Find points to exclude (i.e. outside the bounding polygon)
         from anuga.geometry.polygon import outside_polygon
         cut_points = outside_polygon(gridXY_array, bounding_polygon)
 
+    hole_points_list = []
+    if internal_holes is not None:
+        # Find points to exclude (i.e. inside the internal_holes)
+        from anuga.geometry.polygon import inside_polygon
+        for hole in internal_holes:
+            cut_holes = inside_polygon(gridXY_array, hole)
+            hole_points_list.append(cut_holes)
+
     # Loop over all output quantities and produce the output
     for myTSindex, myTSi in enumerate(myTimeStep):
         if(verbose):
@@ -1150,6 +1160,11 @@ def myInterpFun(quantity):
                 # Cut the points outside the bounding polygon
                 gridq[cut_points] = numpy.nan
 
+            if (internal_holes is not None) and (len(hole_points_list[0]) > 0):
+                # Cut the points inside the hole polygons
+                for hole_points in hole_points_list:
+                    gridq[hole_points] = numpy.nan
+
             # Make name for output file
             if(myTS != 'pointData'):
                 output_name = output_dir + '/' +\