Context Navigation

← Previous Changeset
Next Changeset →

Changeset 9320

Timestamp:

Sep 3, 2014, 3:47:05 PM (11 years ago)

Author:

davies

Message:

Adding a template tsunami script + associated source adjustments

Location:

Files:

: 29 added
: 3 edited

Legend:

: Unmodified
: Added
: Removed

trunk/anuga_core/source/anuga/shallow_water/swDE1_domain_ext.c

r9318	r9320
2092	2092
2093	2093	double timestep;
2094		~~int max_flux_update_frequency;~~
2095	2094
2096	2095	if (!PyArg_ParseTuple(args, "Od", &domain, &timestep)) {

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

-                      r9299
+                      r9320
 """
 import copy
+import os
+import anuga.utilities.spatialInputUtil as su
 def make_nearestNeighbour_quantity_function(quantity_xyValueIn, domain, threshold_distance = 9.0e+100, background_value = 9.0e+100):
 …
         poly_fun_pairs = [ [p0, f0], [p1, f1], ...]
+                    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)
+                    Where:
+                      fi is a function, or a constant, or the name of a
+                        gdal-compatible rasterFile, or a numpy array with 3 columns;
+                      pi is a polygon, or a polygon filename (shapefile or a csv format that anuga.read_polygon will read),
+                        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.
 …
                   pi are polygons where we want to use fi inside
                   SPECIAL pi CASES:
+                  If pi is a filename ending in .shp or a csv format that anuga.read_polygon can read, we assume it contains a polygon we have to read
                   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
 …
                     # Here fi MUST be a gdal-compatible raster
                     # Then we get the extent from the raster itself
-                    import anuga.utilities.spatialInputUtil as su
                     pi_path=su.getRasterExtent(fi,asPolygon=True)
+                elif( type(pi)==str and os.path.isfile(pi) ):
+                    # pi is a file
+                    pi_path=su.read_polygon(pi)
                 else:
+                    # pi is the actual polygon data
                     pi_path=pi
                 notSet=(isSet==0.).nonzero()[0]
                 fInds = inside_polygon(xy_array_trans[notSet,:], pi_path)

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

-                      r9299
+                      r9320
     readShp_1PolyGeo, readShp_1LineGeo -- read SIMPLE shapefile geometries into ANUGA as a list of points.
                                           The supported geometries are quite restricted, see help
+    read_polygon -- Reads either a polygon or line shapefile or a csv as a polygon
     readShpPointsAndAttributes -- read a multi-point shapefile with its attributes into ANUGA
 …
 import sys
 import os
+import os.path
 import copy
 import struct
 import numpy
 #from matplotlib import path
+import anuga
 from anuga.geometry.polygon import inside_polygon
 …
         #dataSrc=ogr.Open(shapefile)
         layer=dataSrc.GetLayer()
+        # Check it is a polygon
+        layerType=ogr.GeometryTypeToName(layer.GetGeomType())
+        if not layerType=='Polygon':
+            msg= shapefile +' is not a polygon shapefile'
+            raise Exception, msg
         # Need a single polygon
         try:
 …
         #dataSrc=ogr.Open(shapefile)
         layer=dataSrc.GetLayer()
+        # Check it is a line
+        layerType=ogr.GeometryTypeToName(layer.GetGeomType())
+        if not layerType=='Line String':
+            msg= shapefile +' is not a line shapefile'
+            raise Exception, msg
         # Need a single line
         try:
             assert(len(layer)==1)
+            assert len(layer)==1
         except:
             print shapefile
 …
     ####################
+    def read_polygon(filename):
+        """
+        Read a shapefile (polygon or line) or an anuga_polygon csv file as an anuga polygon
+        Try to automatically do the correct thing based on the filename
+        """
+        # Check the file exists
+        msg= 'Could not read '+ filename
+        assert os.path.isfile(filename), msg
+        # Get the file extension type
+        fileNameNoExtension , fileExtension = os.path.splitext(filename)
+        if fileExtension == '.shp':
+            # Read as either a polygon or line shapefile
+            try:
+                outPol=readShp_1PolyGeo(filename)
+                assert len(outPol)>1
+            except:
+                try:
+                    outPol= readShp_1LineGeo(filename)
+                    assert len(outPol)>1
+                except:
+                    msg= 'Could not read '+ filename +' as either polygon or line shapefile'
+        else:
+            try:
+                # Read as an anuga polygon file
+                outPol=anuga.read_polygon(filename)
+            except:
+                msg = 'Failed reading polygon '+ filename + ' with anuga.utilities.spatialInputUtils.read_polygon'
+                raise Exception, msg
+        return outPol
+    ####################
     def readShpPtsAndAttributes(shapefile):
         """
 …
         return [pts, attribute, attributeNames]
+    ########################################
+    def readShpPts(shapefile):
+        """
+            Wrapper around readShpPtsAndAttributes
+            Only get the points
+        """
+        out=readShpPtsAndAttributes(shapefile)[0]
+        return out
+    ########################################
+    def read_points(filename):
+        """
+            Reads x,y geometries from a point shapefile or csv file (anuga polygon type format),
+            and returns as a list of lists
+        """
+        # Check the file exists
+        msg= 'Could not read '+ filename
+        assert os.path.isfile(filename), msg
+        # Get the file extension type
+        fileNameNoExtension , fileExtension = os.path.splitext(filename)
+        if fileExtension=='.shp':
+            try:
+                points = readShpPts(filename)
+            except:
+                msg = 'Could not read points from ' + filename
+                raise Exception, msg
+        else:
+            # Assume txt format
+            try:
+                #points=numpy.genfromtxt(filename,delimiter=',',skip_header=1)
+                #points=points[:,0:2].tolist()
+                points=anuga.read_polygon(filename)
+            except:
+                msg = 'Could not read points from ' + filename +\
+                      '. Make sure it has a single header row, with comma separator, and the first 2 columns are x,y'
+                raise Exception, msg
+        return points
     ########################################
 …
         breakLines=copy.copy(breakLinesIn)
         riverWalls=copy.copy(riverWallsIn)
+        # Quick exit
+        if (breakLines == {}) and (riverWalls == {}):
+            return [bounding_polygon, breakLines, riverWalls]
         # Clean intersections of breakLines with itself
         if(verbose):
 …
     #########################################
     def readListOfBreakLines(shapefileList):
         """
             Take a list with the names of shapefiles
+    def readListOfBreakLines(fileList):
+        """
+            Take a list with the names of shapefiles or anuga_polygon csv files
             They are assumed to be '2D breaklines', so we just read their
 …
             Read them in
             INPUT: shapefileList -- a list of shapefile names [e.g. from glob.glob('GIS/Breaklines/*.shp')]
+            INPUT: fileList -- a list of shapefile and/or anuga_polygon csv filenames [e.g. from glob.glob('GIS/Breaklines/*.shp')]
             OUTPUT: dictionary with breaklines [filenames are keys]
 …
         allBreakLines={}
         for shapefile in shapefileList:
             allBreakLines[shapefile]=readShp_1LineGeo(shapefile)
+        for shapefile in fileList:
+            allBreakLines[shapefile]=read_polygon(shapefile) #readShp_1LineGeo(shapefile)
         return allBreakLines

Note: See TracChangeset for help on using the changeset viewer.

Download in other formats: