Changeset 2520

Timestamp:

Mar 10, 2006, 9:07:22 AM (19 years ago)

Author:

sexton

Message:

fixed dem2pts AGAIN!!! to remove NODATA values when a clipping in place

Location:

inundation/pyvolution

Files:

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

inundation/pyvolution/data_manager.py

@@ -1200 +1200 @@
     totalnopoints = nrows*ncols
 
-    # calculating number of NODATA_values for each row
+    # calculating number of NODATA_values for each row in clipped region
+    #FIXME: use array operations to do faster
     nn = 0
+    indices = []
     for i in range(nrows):
         y = (nrows-i)*cellsize + yllcorner
@@ -1209 +1211 @@
                northing_min <= y <= northing_max and \
                dem_elevation_r[i,j] == NODATA_value:
-                nn += 1
+                nn += 1 
+
+            if easting_min <= x <= easting_max and \
+               northing_min <= y <= northing_max:
+                indices.append([i,j])
+
+    # finding indices of the clipped region
+    if len(indices) > 0:
+        firstpoint = indices[0]
+        index1 = firstpoint[0]
+        lastpoint = indices[len(indices)-1]
+        index2 = lastpoint[0]+1 #?
+    else:
+        index1 = 1
+        index2 = ncols
         
     # dimension definitions
@@ -1222 +1238 @@
         print 'There are %d values in the clipped elevation' %clippednopoints
         print 'There are %d NODATA_values in the clipped elevation' %nn
-   
+        
     outfile.createDimension('number_of_points', nopoints)
     outfile.createDimension('number_of_dimensions', 2) #This is 2d data
@@ -1235 +1251 @@
     elevation = outfile.variables['elevation']
 
+    
     #Store data
     global_index = 0
@@ -1242 +1259 @@
 
         lower_index = global_index
-       
-        v = dem_elevation_r[i,:]
+
+        v = dem_elevation_r[i,index1:index2]
         no_NODATA = sum(v == NODATA_value)
         if no_NODATA > 0:
@@ -1249 +1266 @@
         else:
             newcols = ncols_in_bounding_box
-        
+
         telev = zeros(newcols, Float)
         tpoints = zeros((newcols, 2), Float)

inundation/pyvolution/test_data_manager.py

@@ -817 +817 @@
         fid.close()
 
-        nn = ref_elevation.count(NODATA_value)
-        
+       
         #Write prj file with metadata
         metafilename = root+'.prj'
@@ -867 +866 @@
         os.remove(root + '.prj')
 
+    def test_dem2pts_bounding_box_Nullvalues(self):
+        """Test conversion from dem in ascii format to native NetCDF xya format
+        """
+
+        import time, os
+        from Numeric import array, zeros, allclose, Float, concatenate
+        from Scientific.IO.NetCDF import NetCDFFile
+
+        #Write test asc file
+        root = 'demtest'
+
+        filename = root+'.asc'
+        fid = open(filename, 'w')
+        fid.write("""ncols         5
+nrows         6
+xllcorner     2000.5
+yllcorner     3000.5
+cellsize      25
+NODATA_value  -9999
+""")
+        #Create linear function
+
+        ref_points = []
+        ref_elevation = []
+        new_ref_pts1 = []
+        new_ref_elev1 = []
+        NODATA_value = -9999
+        for i in range(6):
+            y = (6-i)*25.0
+            for j in range(5):
+                x = j*25.0
+                z = x+2*y
+                if j == 4: z = NODATA_value # column
+                if i == 2 and j == 2: z = NODATA_value # random
+                if i == 5 and j == 1: z = NODATA_value
+                if i == 1: z = NODATA_value # row
+                if i == 3 and j == 1: z = NODATA_value # two pts/row
+                if i == 3 and j == 3: z = NODATA_value
+
+                if z <> NODATA_value:
+                    new_ref_elev1.append(z)
+                    new_ref_pts1.append( [x,y] )
+                    
+                ref_points.append( [x,y] )
+                ref_elevation.append(z)
+                fid.write('%f ' %z)
+            fid.write('\n')
+
+        fid.close()
+
+        #Write prj file with metadata
+        metafilename = root+'.prj'
+        fid = open(metafilename, 'w')
+
+
+        fid.write("""Projection UTM
+Zone 56
+Datum WGS84
+Zunits NO
+Units METERS
+Spheroid WGS84
+Xshift 0.0000000000
+Yshift 10000000.0000000000
+Parameters
+""")
+        fid.close()
+       
+        #Convert to NetCDF pts
+        convert_dem_from_ascii2netcdf(root)
+        dem2pts(root, easting_min=2010.0, easting_max=2110.0,
+                northing_min=3035.0, northing_max=3125.5)
+
+        #Check contents
+        #Get NetCDF
+        fid = NetCDFFile(root+'.pts', 'r')
+
+        # Get the variables
+        #print fid.variables.keys()
+        points = fid.variables['points']
+        elevation = fid.variables['elevation']
+
+        #Check values
+        assert fid.xllcorner[0] == 2010.0
+        assert fid.yllcorner[0] == 3035.0
+
+        #create new reference points
+        ref_points = []
+        ref_elevation = []
+        new_ref_pts2 = []
+        new_ref_elev2 = []
+        for i in range(4):
+            y = (4-i)*25.0 + 25.0
+            y_new = y + 3000.5 - 3035.0
+            for j in range(4):
+                x = j*25.0 + 25.0
+                x_new = x + 2000.5 - 2010.0
+                z = x+2*y
+
+                if j == 3: z = NODATA_value # column
+                if i == 1 and j == 1: z = NODATA_value # random
+                if i == 4 and j == 0: z = NODATA_value
+                if i == 0: z = NODATA_value # row
+                if i == 2 and j == 0: z = NODATA_value # two pts/row
+                if i == 2 and j == 2: z = NODATA_value
+
+                if z <> NODATA_value:
+                    new_ref_elev2.append(z)
+                    new_ref_pts2.append( [x_new,y_new] )
+
+                    
+                ref_points.append( [x_new,y_new] )
+                ref_elevation.append(z)
+
+        #print points[:]
+        #print ref_points
+        #assert allclose(points, ref_points)
+
+        #print attributes[:]
+        #print ref_elevation
+        #assert allclose(elevation, ref_elevation)
+
+        
+        assert len(points) == len(new_ref_pts2), 'length of returned points not correct'
+        assert allclose(points, new_ref_pts2), 'points do not align'
+
+        #print 'elevation', elevation[:]
+        assert len(elevation) == len(new_ref_elev2), 'length of returned elevation not correct'
+        assert allclose(elevation, new_ref_elev2), 'elevations do not align'
+        #Cleanup
+        fid.close()
+
+
+        os.remove(root + '.pts')
+        os.remove(root + '.dem')
+        os.remove(root + '.asc')
+        os.remove(root + '.prj')