Changeset 1875

Timestamp:

Oct 6, 2005, 2:38:17 PM (19 years ago)

Author:

ole

Message:

Bounding box for sww2dem

Location:

inundation/pyvolution

Files:

• data_manager.py (modified) (9 diffs)
• test_data_manager.py (modified) (6 diffs)

inundation/pyvolution/data_manager.py

@@ -1211 +1211 @@
             reduction = None,
             cellsize = 10,
+            easting_min = None,
+            easting_max = None,
+            northing_min = None,
+            northing_max = None,          
             verbose = False,
             origin = None,
@@ -1277 +1281 @@
     # Note the use of a .ers extension is optional (write_ermapper_grid will
     # deal with either option
-    # ersfile = basename_out
-
+
+    #Read sww file
     if verbose: print 'Reading from %s' %swwfile
-    #Read sww file
     from Scientific.IO.NetCDF import NetCDFFile
     fid = NetCDFFile(swwfile)
@@ -1288 +1291 @@
     y = fid.variables['y'][:]
     volumes = fid.variables['volumes'][:]
-
-    ymin = min(y); ymax = max(y)
-    xmin = min(x); xmax = max(x)
 
     number_of_timesteps = fid.dimensions['number_of_timesteps']
@@ -1312 +1312 @@
 
 
+
     #Get quantity and reduce if applicable
     if verbose: print 'Reading quantity %s' %quantity
@@ -1322 +1323 @@
 
     if len(q.shape) == 2:
+        #q has a time component and needs to be reduced along
+        #the temporal dimension
         if verbose: print 'Reducing quantity %s' %quantity
         q_reduced = zeros( number_of_points, Float )
@@ -1330 +1333 @@
         q = q_reduced
 
-    #Now q has dimension: number_of_points
+    #Post condition: Now q has dimension: number_of_points
+    assert len(q.shape) == 1
+    assert q.shape[0] == number_of_points    
+    
 
     #Create grid and update xll/yll corner and x,y
+
+    #Relative extent
+    if easting_min is None:
+        xmin = min(x)
+    else:
+        xmin = easting_min - xllcorner
+
+    if easting_max is None:
+        xmax = max(x)
+    else:
+        xmax = easting_max - xllcorner
+        
+    if northing_min is None:        
+        ymin = min(y)
+    else:
+        ymin = northing_min - yllcorner
+        
+    if northing_max is None:        
+        ymax = max(y)
+    else:
+        ymax = northing_max - yllcorner
+        
+        
+    
     if verbose: print 'Creating grid'
     ncols = int((xmax-xmin)/cellsize)+1
     nrows = int((ymax-ymin)/cellsize)+1
 
+
+    #New absolute reference and coordinates
     newxllcorner = xmin+xllcorner
     newyllcorner = ymin+yllcorner
@@ -1371 +1403 @@
     #FIXME: This should be done with precrop = True (?), otherwise it'll
     #take forever. With expand_search  set to False, some grid points might
-    #miss out....
+    #miss out.... This will be addressed though Duncan's refactoring of least_squares
 
     interp = Interpolation(vertex_points, volumes, grid_points, alpha=0.0,
@@ -1474 +1506 @@
             origin = None):
     print 'sww2asc will soon be obsoleted - please use sww2dem'
-    sww2dem(basename_in, 
-            basename_out,
-            quantity,
-            timestep,
-            reduction,
-            cellsize,
-            verbose,
-            origin,
+    sww2dem(basename_in,
+            basename_out = basename_out,
+            quantity = quantity,
+            timestep = timestep,
+            reduction = reduction,
+            cellsize = cellsize,
+            verbose = verbose,
+            origin = origin,
             datum = 'WGS84',
             format = 'asc')
@@ -1495 +1527 @@
     print 'sww2ers will soon be obsoleted - please use sww2dem'
     sww2dem(basename_in, 
-            basename_out,
-            quantity,
-            timestep,
-            reduction,
-            cellsize,
-            verbose,
-            origin,
-            datum,
+            basename_out = basename_out,
+            quantity = quantity,
+            timestep = timestep,
+            reduction = reduction,
+            cellsize = cellsize,
+            verbose = verbose,
+            origin = origin,
+            datum = datum,
             format = 'ers')         
 ################################# END COMPATIBILITY ##############          

inundation/pyvolution/test_data_manager.py

@@ -965 +965 @@
                 verbose = False,
                 format = 'asc')
-        #dem2pts(root, easting_min=2010.0, easting_max=2110.0,
-        #        northing_min=3035.0, northing_max=3125.5)
-        
         
                 
@@ -1056 +1053 @@
 
 
-    #FIXME: TODO
-    def xtest_sww2dem_boundingbox(self):
+    def test_sww2dem_boundingbox(self):
         """Test that sww information can be converted correctly to asc/prj
         format readable by e.g. ArcView.
@@ -1158 +1154 @@
                 format = 'asc')
         
-                
+        fid.close()
+
+        
         #Check prj (meta data)
         prjid = open(prjfile)
@@ -1236 +1234 @@
         
 
-        fid.close()
 
         #Cleanup
@@ -1593 +1590 @@
         os.remove(sww.filename)
 
-
-    def xxxtestz_sww2ers_real(self):
-        """Test that sww information can be converted correctly to asc/prj
-        format readable by e.g. ArcView
-        """
-
-        import time, os
-        from Numeric import array, zeros, allclose, Float, concatenate
-        from Scientific.IO.NetCDF import NetCDFFile
-
-
-
-
-        #Export to ascii/prj files
-        sww2ers('karratha_100m.sww',
-                quantity = 'depth',
-                cellsize = 5,
-                verbose = False)
 
 
@@ -2541 +2520 @@
 if __name__ == "__main__":
     suite = unittest.makeSuite(Test_Data_Manager,'test')
+    #suite = unittest.makeSuite(Test_Data_Manager,'test_sww2dem_boundingbox')    
     #suite = unittest.makeSuite(Test_Data_Manager,'test_dem2pts_bounding_box')
     #suite = unittest.makeSuite(Test_Data_Manager,'test_decimate_dem')