Changeset 9219

Timestamp:

Jun 25, 2014, 4:55:20 PM (10 years ago)

Author:

davies

Message:

Changes + enhanced testing for spatialInputUtil + quantity_setting_functions

Location:

trunk/anuga_core/source/anuga/utilities

Files:

• quantity_setting_functions.py (modified) (8 diffs)
• spatialInputUtil.py (modified) (2 diffs)
• test_quantity_setting_functions.py (modified) (2 diffs)
• test_spatialInputUtil.py (modified) (2 diffs)

trunk/anuga_core/source/anuga/utilities/quantity_setting_functions.py

@@ -78 +78 @@
         Make a 'composite function' to set quantities -- applies different functions inside different polygon regions.
              
-        fun_poly_pairs = [ [f0, p0], [f1, p1], ...] where fi is a function and pi is a polygon, or None, or 'All'
+        fun_poly_pairs = [ [f0, p0], [f1, p1], ...] 
+                    where fi is a function, or a constant, or the name of a gdal-compatible rasterFile; 
+                    and pi is a polygon, or None, or 'All' (or it can be 'Extent' in the case that fi is a rasterFile name)
               
         IMPORTANT: When polygons overlap, the first elements of the list are given priority. 
@@ -90 +92 @@
               fun_poly_pairs = [ [f0, p0], [f1, p1], ...]
 
-                  where fi(x,y) are functions returning quantity values at points, and 
-                  pi are polygons where we want to use fi inside
+                  where fi(x,y) are functions returning quantity values at
+                  points, or constants in which case points in the polygon are set to that
+                  value, or (string) rasterFile names which can be passed to quantityRasterFun to make a function
+
+                  and pi are polygons where we want to use fi inside
                 
                   If any pi = 'All', then we assume that ALL unset points are set
-                    using the function. This CAN ONLY happen in the last [fi,pi] pair where pi is
-                    not None (since fi will be applied to all remaining points -- so anything else is probably an input mistake)
+                  using the function. This CAN ONLY happen in the last [fi,pi] pair where pi is
+                  not None (since fi will be applied to all remaining points -- so anything else is probably an input mistake)
     
                   If any pi = None, then that [fi,pi] pair is skipped
 
+                  if pi='Extent' and fi is the name of a raster file, then the
+                    extent of the raster file is used to define the polygon
+
               domain = ANUGA domain object
 
@@ -105 +113 @@
                 domain.set_quantity('elevation', F)
     """
+    import os
     import scipy
     from matplotlib import path
@@ -136 +145 @@
             fi = fun_poly_pairs[i][0] # The function
             pi = fun_poly_pairs[i][1] # The polygon
+
+            # Quick exit
+            if(pi is None):
+                continue
+
             # Get a 'true/false' flag for whether the point is in pi, and not set
             if(pi is 'All'):
                 fInside=(1-isSet)
-            elif(pi is None):
-                continue
             else:
-                # Try matplotlib -- make the path a matplotlib path object
-                pi_path=path.Path(scipy.array(pi))
+                if(pi is 'Extent' and type(fi) is str and os.path.exists(fi)):
+                    # Here fi MUST be a gdal-compatible raster
+                    # Then we get the extent from the raster itself
+                    import anuga.utilities.spatialInputUtil as su
+                    newpi=su.getRasterExtent(fi,asPolygon=True)
+                    pi_path=path.Path(scipy.array(newpi))
+                else:
+                    # Try matplotlib -- make the path a matplotlib path object
+                    pi_path=path.Path(scipy.array(pi))
+
                 fInside=xy_array_trans[:,0]*0.
                 for j in range(len(fInside)):
@@ -153 +173 @@
             fInds=(fInside==1.).nonzero()[0]
             if(len(fInds)>0):
-                quantityVal[fInds] = fi(x[fInds], y[fInds])
+                if(hasattr(fi,'__call__')):
+                    # fi is a function
+                    quantityVal[fInds] = fi(x[fInds], y[fInds])
+                elif isinstance(fi, (int, long, float)):
+                    # fi is a numerical constant
+                    quantityVal[fInds]=fi*1.0
+                elif ( type(fi) is str and os.path.exists(fi)):
+                    # fi is a file which is assumed to be 
+                    # a gdal-compatible raster
+                    newfi = quantityRasterFun(domain, fi)
+                    quantityVal[fInds] = newfi(x[fInds], y[fInds])
+                
                 isSet[fInds]=1
 
@@ -191 +222 @@
 #################################################################################
 
-def elevation_from_Pt_Pol_Data_and_Raster(Pt_Pol_Data, elevation_raster, domain):
-    """
-        Function to make a function F(x,y) which returns the elevation on elevation_raster,
-            except if x,y is inside the polygon associated with any element of Pt_Pol_Data,
-            in which case a Pt_Pol_-specific nearest neighbour interpolator is used.
+def quantity_from_Pt_Pol_Data_and_Raster(Pt_Pol_Data, quantity_raster, domain):
+    """
+        Function to make a function F(x,y) which returns the corresponding
+        values on quantity_raster, except if x,y is inside the polygon associated with
+        any element of Pt_Pol_Data, in which case a Pt_Pol_-specific nearest neighbour
+        interpolator is used.
 
         INPUT:
@@ -201 +233 @@
                                         [ Polygon_1, Pt_XYZ_1], 
                                         ... ]
-            elevation_raster = A GDAL-compatible elevation raster
+                    Here Polygon_i is a polygon in ANUGA format,
+                    and Pt_XYZ_i is a 3 column array of x,y,Value points
+            quantity_raster = A GDAL-compatible quantity raster
             domain = ANUGA domain
     """
 
-    # Function to set elevation from raster 
-    elevFun1=quantityRasterFun(domain, rasterFile=elevation_raster)
-
-    # List of function/polygon pairs defining the Pt_Pol_ elevation data
-    elevFunChanList=[]
+    # Function to set quantity from raster 
+    qFun1=quantityRasterFun(domain, rasterFile=quantity_raster)
+
+    # List of function/polygon pairs defining the Pt_Pol_ quantity data
+    qFunChanList=[]
     for i in range(len(Pt_Pol_Data)):
-        elevFunChanList.append([
+        qFunChanList.append([
              make_nearestNeighbour_quantity_function(Pt_Pol_Data[i][1], domain), 
              Pt_Pol_Data[i][0]
@@ -217 +251 @@
 
     #
-    elevFun=composite_quantity_setting_function(elevFunChanList+[[elevFun1, 'All']], domain)
-
-    return elevFun
+    qFun=composite_quantity_setting_function(qFunChanList+[[qFun1, 'All']], domain)
+
+    return qFun

trunk/anuga_core/source/anuga/utilities/spatialInputUtil.py

@@ -502 +502 @@
         else:
             return [L1, L2]
-    
+   
+    ###########################################################
+    def getRasterExtent(rasterFile, asPolygon=False):
+        """
+            Sometimes we need to know the extent of a raster
+            i.e. the minimum x, maximum x, minimum y, and maximum y values
+            
+            INPUT:
+                rasterFile -- a gdal compatible rasterfile
+                asPolygon -- if False, return [xmin,xmax,ymin,ymax].
+                             If True, return [ [xmin,ymin],[xmax,ymin],[xmax,ymax],[xmin,ymax]]
+            OUTPUT
+                The extent as defined above
+
+        """
+        raster = gdal.Open(rasterFile)
+        transform=raster.GetGeoTransform()
+        xOrigin = transform[0]
+        yOrigin = transform[3]
+        xPixels=raster.RasterXSize
+        yPixels=raster.RasterYSize
+
+        # Compute the other extreme corner
+        x2=xOrigin + xPixels*transform[1]+yPixels*transform[2]
+        y2=yOrigin + xPixels*transform[4]+yPixels*transform[5]
+        
+        xmin=min(xOrigin,x2) 
+        xmax=max(xOrigin,x2)
+
+        ymin=min(yOrigin,y2)
+        ymax=max(yOrigin,y2)
+
+        if(asPolygon):
+            return [ [xmin,ymin], [xmax,ymin], [xmax,ymax], [xmin,ymax]]
+        else:
+            return [xmin,xmax,ymin,ymax]
+
+ 
     ###########################################################
     def rasterValuesAtPoints(xy, rasterFile, band=1):
@@ -548 +585 @@
             print 'You might need to edit this code to read the data type'
             raise Exception, 'Stopping'
-    
+  
+        # Upper bounds for pixel values, so we can fail gracefully
+        xMax=raster.RasterXSize
+        yMax=raster.RasterYSize
+        if(px.max()<xMax and px.min()>=0 and py.max()<yMax and py.min()>=0):
+            pass
+        else:
+            raise Exception, 'Trying to extract point values that exceed the raster extent'
+
         # Get values -- seems we have to loop, but it is efficient enough
         for i in range(len(px)):

trunk/anuga_core/source/anuga/utilities/test_quantity_setting_functions.py

@@ -140 +140 @@
         return
 
-    def test_elevation_from_Pt_Pol_Data_and_Raster(self):
-        # 
-        # 
-        #
+    def test_composite_quantity_setting_function(self):
+        # Test the composite_quantity_setting_function
+        
         domain=self.create_domain(1.0, 0.0)
-
-        # Evolve the model
-        #for t in domain.evolve(yieldstep=0.2, finaltime=1.0):
-        #    pass
-
+        
         # Make a raster from the elevation data
         from anuga.utilities import plot_utils as util
@@ -162 +157 @@
         # Make a polygon-point pair which we use to set elevation in a 'channel'
         trenchPoly=[[minX+40., minY], [minX+40., minY+100.], [minX+60., minY+100.], [minX+60., minY]]
+
+        #################################################################
+ 
+        # This example uses a constant, and a raster, to set the quantity           
+        F=qs.composite_quantity_setting_function([[-1000., trenchPoly], ['PointData_ElevTest.tif', 'Extent']],\
+                                                domain) 
+
+        # Points where we test the function
+        testPts_X=numpy.array([50., 3.])
+        testPts_Y=numpy.array([1., 20.])
+        fitted=F(testPts_X,testPts_Y)
+
+        # The fitted value in the trench should be -1000.
+        assert(fitted[0]==-1000.)
+
+        # Find the nearest domain point to the second test point
+        # This will have been used in constructing the elevation raster
+        nearest=((domain.centroid_coordinates[:,0]-3.)**2 + (domain.centroid_coordinates[:,1]-20.)**2).argmin()
+        nearest_x=domain.centroid_coordinates[nearest,0]
+        assert(numpy.allclose(fitted[1],-nearest_x/150.))
+
+        #########################################################################
+
+        # This example uses a function, and a raster, to set the quantity           
+        def f0(x,y):
+            return x/10.
+        F=qs.composite_quantity_setting_function([[f0, trenchPoly], ['PointData_ElevTest.tif', 'Extent']],\
+                                                domain) 
+        fitted=F(testPts_X,testPts_Y)
+        # Now the fitted value in the trench should be determined by f0
+        assert(numpy.allclose(fitted[0],50./10.))
+        # The second test point should be as before
+        nearest=((domain.centroid_coordinates[:,0]-3.)**2 + (domain.centroid_coordinates[:,1]-20.)**2).argmin()
+        nearest_x=domain.centroid_coordinates[nearest,0]
+        assert(numpy.allclose(fitted[1],-nearest_x/150.))
+
+        ##########################################################################
+
+        # This example uses 'All' as a polygon
+        F=qs.composite_quantity_setting_function([[f0, 'All'], ['PointData_ElevTest.tif', None]],\
+                                                domain) 
+        fitted=F(testPts_X,testPts_Y)
+        # Now the fitted value in the trench should be determined by f0
+        assert(numpy.allclose(fitted[0],50./10.))
+        assert(numpy.allclose(fitted[1],3./10.))
+
+        ###########################################################################
+        # This example should fail
+        try:
+            F=qs.composite_quantity_setting_function([[f0, 'All'], ['PointData_ElevTest.tif', 'All']],\
+                                                    domain) 
+            raise Exception, 'The last command should fail' 
+        except:
+            assert True
+
+        return
+
+    def test_quantity_from_Pt_Pol_Data_and_Raster(self):
+        # 
+        # 
+        #
+        domain=self.create_domain(1.0, 0.0)
+
+        # Evolve the model
+        #for t in domain.evolve(yieldstep=0.2, finaltime=1.0):
+        #    pass
+
+        # Make a raster from the elevation data
+        from anuga.utilities import plot_utils as util
+        xs=domain.centroid_coordinates[:,0]+domain.geo_reference.xllcorner
+        ys=domain.centroid_coordinates[:,1]+domain.geo_reference.yllcorner
+        elev=domain.quantities['elevation'].centroid_values
+
+        allDat=numpy.vstack([xs,ys,elev]).transpose()
+        util.Make_Geotif(allDat, output_quantities=['ElevTest'], EPSG_CODE=32756, 
+                        output_dir='.', CellSize=1.)
+
+        # Make a polygon-point pair which we use to set elevation in a 'channel'
+        trenchPoly=[[minX+40., minY], [minX+40., minY+100.], [minX+60., minY+100.], [minX+60., minY]]
         trenchPts=numpy.array([minX+50., minY+50., -1000.])
         #
         PtPolData=[[trenchPoly, trenchPts]]
-        F=qs.elevation_from_Pt_Pol_Data_and_Raster(PtPolData, 'PointData_ElevTest.tif', domain) 
+        F=qs.quantity_from_Pt_Pol_Data_and_Raster(PtPolData, 'PointData_ElevTest.tif', domain) 
 
         testPts_X=numpy.array([50., 3.])

trunk/anuga_core/source/anuga/utilities/test_spatialInputUtil.py

@@ -29 +29 @@
     def tearDown(self):
         pass
+
+    def make_me_a_tif(self):
+        # We need to make a .tif to test some functions 
+        # This does the job
+        #
+        from anuga.utilities import plot_utils as util
+        #
+        # Do it with Make_Geotif
+        # Pick a domain that makes sense in EPSG:32756
+        x=numpy.linspace(307000., 307100., 101)     
+        y=numpy.linspace(6193000., 6193100., 101)
+        xG,yG=numpy.meshgrid(x, y)
+        xG=xG.flatten()
+        yG=yG.flatten()
+        # Surface is z=x+y
+        fakeZ=xG-min(xG)+yG -min(yG)
+        dataToGrid=numpy.vstack([xG,yG,fakeZ]).transpose()
+        #    
+        util.Make_Geotif(dataToGrid, output_quantities=['TestData'],
+                         EPSG_CODE=32756, output_dir='.',CellSize=1.0)
    
     def test_compute_square_distance_to_segment(self):
@@ -281 +301 @@
         assert(all(pip[:,0]+pip[:,1]>=10.))
 
+
+    def test_getRasterExtent(self):
+        self.make_me_a_tif()
+
+        extentOut=su.getRasterExtent('PointData_TestData.tif')
+        assert(numpy.allclose(extentOut[0], 307000.-0.5))
+        assert(numpy.allclose(extentOut[1], 307100.+0.5))
+        assert(numpy.allclose(extentOut[2], 6193000.-0.5))
+        assert(numpy.allclose(extentOut[3], 6193100.+0.5))
+        
+        extentOut=su.getRasterExtent('PointData_TestData.tif',asPolygon=True)
+        assert(numpy.allclose(extentOut[0][0], 307000.-0.5))
+        assert(numpy.allclose(extentOut[3][0], 307000.-0.5))
+        assert(numpy.allclose(extentOut[1][0], 307100.+0.5))
+        assert(numpy.allclose(extentOut[2][0], 307100.+0.5))
+        assert(numpy.allclose(extentOut[0][1], 6193000.-0.5))
+        assert(numpy.allclose(extentOut[1][1], 6193000.-0.5))
+        assert(numpy.allclose(extentOut[2][1], 6193100.+0.5))
+        assert(numpy.allclose(extentOut[3][1], 6193100.+0.5))
+
     def test_rasterValuesAtPoints(self):
         # We need to make a .tif to test this function. 
-        from anuga.utilities import plot_utils as util
-        # Do it with Make_Geotif
-        # Pick a domain that makes sense in EPSG:32756
-        x=numpy.linspace(307000., 307100., 101)     
-        y=numpy.linspace(6193000., 6193100., 101)
-        xG,yG=numpy.meshgrid(x, y)
-        xG=xG.flatten()
-        yG=yG.flatten()
-        # Surface is z=x+y
-        fakeZ=xG-min(xG)+yG -min(yG)
-        dataToGrid=numpy.vstack([xG,yG,fakeZ]).transpose()
-
-        util.Make_Geotif(dataToGrid, output_quantities=['TestData'],
-                         EPSG_CODE=32756, output_dir='.',CellSize=1.0)
-
+        self.make_me_a_tif()
+
+        # Get the range of the tif
+        tifRange=su.getRasterExtent('PointData_TestData.tif')
 
         # Now try to get some point values -- note they will be rounded to the
         # nearest cell
-        xA=numpy.array([0., 10.3, 50.9, 100.])+x.min()
-        yA=numpy.array([0., 20.1, 75.1, 100.])+y.min()
-        z_predicted=numpy.round(xA)+numpy.round(yA)-x.min()-y.min()
+        xA=numpy.array([0., 10.3, 50.9, 100.])+tifRange[0]+0.5
+        yA=numpy.array([0., 20.1, 75.1, 100.])+tifRange[2]+0.5
+        z_predicted=numpy.round(xA)+numpy.round(yA)-tifRange[0]-tifRange[2]-1.0
         InDat=numpy.vstack([xA,yA]).transpose()
         z_fitted=su.rasterValuesAtPoints(InDat, rasterFile='PointData_TestData.tif')