Changeset 1577 for inundation/ga/storm_surge/pyvolution/data_manager.py

Timestamp:

Jul 5, 2005, 12:16:43 PM (19 years ago)

Author:

hitchman

Message:

added bounding box to dem2pts + test

File:

• inundation/ga/storm_surge/pyvolution/data_manager.py (modified) (4 diffs)

inundation/ga/storm_surge/pyvolution/data_manager.py

@@ -921 +921 @@
     return cls(domain, mode)
 
-def dem2pts(basename_in, basename_out=None, verbose=False):
+def dem2pts(basename_in, basename_out=None, verbose=False,
+            easting_min=None, easting_max=None,
+            northing_min=None, northing_max=None):
     """Read Digitial Elevation model from the following NetCDF format (.dem)
 
@@ -983 +985 @@
     outfile.description = 'NetCDF pts format for compact and portable storage ' +\
                       'of spatial point data'
+    #assign default values
+    if easting_min is None: easting_min = xllcorner
+    if easting_max is None: easting_max = xllcorner + ncols*cellsize
+    if northing_min is None: northing_min = yllcorner
+    if northing_max is None: northing_max = yllcorner + nrows*cellsize
+
+    #compute offsets to update georeferencing
+    easting_offset = xllcorner - easting_min
+    northing_offset = yllcorner - northing_min
 
     #Georeferencing
     outfile.zone = zone
-    outfile.xllcorner = xllcorner #Easting of lower left corner
-    outfile.yllcorner = yllcorner #Northing of lower left corner
+    outfile.xllcorner = easting_min #Easting of lower left corner
+    outfile.yllcorner = northing_min #Northing of lower left corner
     outfile.false_easting = false_easting
     outfile.false_northing = false_northing
@@ -1002 +1013 @@
 
     # dimension definitions
-    outfile.createDimension('number_of_points', nrows*ncols)
+    nrows_in_bounding_box = int(round((easting_max-easting_min)/cellsize))
+    ncols_in_bounding_box = int(round((northing_max-northing_min)/cellsize))
+    outfile.createDimension('number_of_points', nrows_in_bounding_box*ncols_in_bounding_box)
     outfile.createDimension('number_of_dimensions', 2) #This is 2d data
 
@@ -1016 +1029 @@
     #Store data
     #FIXME: Could perhaps be faster using array operations
+    index = 0
     for i in range(nrows):
         if verbose: print 'Processing row %d of %d' %(i, nrows)
 
-        y = (nrows-i)*cellsize
+        y = (nrows-i)*cellsize + yllcorner
         for j in range(ncols):
-            index = i*ncols + j
-
-            x = j*cellsize
-            points[index, :] = [x,y]
-            elevation[index] = dem_elevation[i, j]
+            #index = i*ncols + j
+
+            x = j*cellsize + xllcorner
+            if easting_min <= x <= easting_max and \
+               northing_min <= y <= northing_max:
+                points[index, :] = [x-easting_min,y-northing_min]
+                elevation[index] = dem_elevation[i, j]
+                index += 1
+
+    #print easting_min, easting_max
+    #print index, nrows_in_bounding_box*ncols_in_bounding_box
+    #print nrows, ncols
+    assert index == nrows_in_bounding_box*ncols_in_bounding_box, 'index not equal to number of points'
 
     infile.close()