Changeset 2060

Timestamp:

Nov 24, 2005, 10:46:22 AM (19 years ago)

Author:

ole

Message:

Made asc and ers formats use the same method for assigning missing values

Location:

inundation/pyvolution

Files:

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

inundation/pyvolution/data_manager.py

@@ -1424 +1424 @@
             reduction = None,
             cellsize = 10,
+            NODATA_value = -9999,
             easting_min = None,
             easting_max = None,
@@ -1481 +1482 @@
     msg = 'Format must be either asc or ers'
     assert format.lower() in ['asc', 'ers'], msg
-    
+
+        
     false_easting = 500000
     false_northing = 10000000
@@ -1670 +1672 @@
     grid_values = interp.interpolate(q).flat
 
+
+    #Assign NODATA_value to all points outside bounding polygon (from interpolation mesh)
+    P = interp.mesh.get_boundary_polygon()
+    inside_indices = inside_polygon(grid_points, P)
+
+    #FIXME: Remove when new interpolation takes care of this
+    for i in range(nrows):
+        if verbose and i%((nrows+10)/10)==0:
+            print 'Doing row %d of %d' %(i, nrows)
+
+        base_index = (nrows-i-1)*ncols     
+        for j in range(ncols):
+            index = base_index+j
+
+            if not sometrue(inside_indices == index):
+                grid_values[index] = NODATA_value
+
+    #############            
+
+    
+
     if verbose:
         print 'Interpolated values are in [%f, %f]' %(min(grid_values),
@@ -1677 +1700 @@
         # setup ERS header information
         grid_values = reshape(grid_values,(nrows, ncols))    
-        NODATA_value = 0
         header = {}
         header['datum'] = '"' + datum + '"'
@@ -1725 +1747 @@
     
         if verbose: print 'Writing %s' %demfile
-        NODATA_value = -9999
  
         ascid = open(demfile, 'w')
@@ -1738 +1759 @@
 
         #Get bounding polygon from mesh
-        P = interp.mesh.get_boundary_polygon()
-        inside_indices = inside_polygon(grid_points, P)
+        #P = interp.mesh.get_boundary_polygon()
+        #inside_indices = inside_polygon(grid_points, P)
 
         for i in range(nrows):
@@ -1745 +1766 @@
                 print 'Doing row %d of %d' %(i, nrows)
 
+            base_index = (nrows-i-1)*ncols     
             for j in range(ncols):
-                index = (nrows-i-1)*ncols+j
-
-                if sometrue(inside_indices == index):
-                    ascid.write('%f ' %grid_values[index])
-                else:
-                    ascid.write('%d ' %NODATA_value)
+                index = base_index+j
+
+                ascid.write('%f ' %grid_values[index])
 
             ascid.write('\n')

inundation/pyvolution/test_data_manager.py

@@ -1694 +1694 @@
         from Scientific.IO.NetCDF import NetCDFFile
 
+
+        NODATA_value = 1758323
+
         #Setup
         self.domain.filename = 'datatest'
@@ -1737 +1740 @@
                 quantity = 'elevation',
                 cellsize = cellsize,
+                NODATA_value = NODATA_value,
                 verbose = False,
                 format = 'ers')
@@ -1755 +1759 @@
         assert int(header['nroflines']) == 5
         assert int(header['nrofcellsperline']) == 5     
-        assert int(header['nullcellvalue']) == 0   #?           
+        assert int(header['nullcellvalue']) == NODATA_value
         #FIXME - there is more in the header                    
 
@@ -1763 +1767 @@
         
         #FIXME (Ole): Why is this the desired reference grid for -x-y?
-        ref_grid = [0,      0,     0,     0,     0,
+        ref_grid = [NODATA_value, NODATA_value, NODATA_value, NODATA_value, NODATA_value,
                     -1,    -1.25, -1.5,  -1.75, -2.0,
                     -0.75, -1.0,  -1.25, -1.5,  -1.75,              
@@ -1770 +1774 @@
                                           
 
+        #print grid
         assert allclose(grid, ref_grid)
 
@@ -3581 +3586 @@
 if __name__ == "__main__":
     #suite = unittest.makeSuite(Test_Data_Manager,'test_asc_csiro2sww4')
-    suite = unittest.makeSuite(Test_Data_Manager,'test')
-    #suite = unittest.makeSuite(Test_Data_Manager,'xxxtest')
+    #suite = unittest.makeSuite(Test_Data_Manager,'test')
+    suite = unittest.makeSuite(Test_Data_Manager,'test_sww2dem_asc_missing_points')
     #suite = unittest.makeSuite(Test_Data_Manager,'test_sww2dem_boundingbox')    
     #suite = unittest.makeSuite(Test_Data_Manager,'test_dem2pts_bounding_box')