Changeset 2520


Ignore:
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:
2 edited

Legend:

Unmodified
Added
Removed
  • inundation/pyvolution/data_manager.py

    r2518 r2520  
    12001200    totalnopoints = nrows*ncols
    12011201
    1202     # calculating number of NODATA_values for each row
     1202    # calculating number of NODATA_values for each row in clipped region
     1203    #FIXME: use array operations to do faster
    12031204    nn = 0
     1205    indices = []
    12041206    for i in range(nrows):
    12051207        y = (nrows-i)*cellsize + yllcorner
     
    12091211               northing_min <= y <= northing_max and \
    12101212               dem_elevation_r[i,j] == NODATA_value:
    1211                 nn += 1
     1213                nn += 1
     1214
     1215            if easting_min <= x <= easting_max and \
     1216               northing_min <= y <= northing_max:
     1217                indices.append([i,j])
     1218
     1219    # finding indices of the clipped region
     1220    if len(indices) > 0:
     1221        firstpoint = indices[0]
     1222        index1 = firstpoint[0]
     1223        lastpoint = indices[len(indices)-1]
     1224        index2 = lastpoint[0]+1 #?
     1225    else:
     1226        index1 = 1
     1227        index2 = ncols
    12121228       
    12131229    # dimension definitions
     
    12221238        print 'There are %d values in the clipped elevation' %clippednopoints
    12231239        print 'There are %d NODATA_values in the clipped elevation' %nn
    1224    
     1240        
    12251241    outfile.createDimension('number_of_points', nopoints)
    12261242    outfile.createDimension('number_of_dimensions', 2) #This is 2d data
     
    12351251    elevation = outfile.variables['elevation']
    12361252
     1253   
    12371254    #Store data
    12381255    global_index = 0
     
    12421259
    12431260        lower_index = global_index
    1244        
    1245         v = dem_elevation_r[i,:]
     1261
     1262        v = dem_elevation_r[i,index1:index2]
    12461263        no_NODATA = sum(v == NODATA_value)
    12471264        if no_NODATA > 0:
     
    12491266        else:
    12501267            newcols = ncols_in_bounding_box
    1251        
     1268
    12521269        telev = zeros(newcols, Float)
    12531270        tpoints = zeros((newcols, 2), Float)
  • inundation/pyvolution/test_data_manager.py

    r2514 r2520  
    817817        fid.close()
    818818
    819         nn = ref_elevation.count(NODATA_value)
    820        
     819       
    821820        #Write prj file with metadata
    822821        metafilename = root+'.prj'
     
    867866        os.remove(root + '.prj')
    868867
     868    def test_dem2pts_bounding_box_Nullvalues(self):
     869        """Test conversion from dem in ascii format to native NetCDF xya format
     870        """
     871
     872        import time, os
     873        from Numeric import array, zeros, allclose, Float, concatenate
     874        from Scientific.IO.NetCDF import NetCDFFile
     875
     876        #Write test asc file
     877        root = 'demtest'
     878
     879        filename = root+'.asc'
     880        fid = open(filename, 'w')
     881        fid.write("""ncols         5
     882nrows         6
     883xllcorner     2000.5
     884yllcorner     3000.5
     885cellsize      25
     886NODATA_value  -9999
     887""")
     888        #Create linear function
     889
     890        ref_points = []
     891        ref_elevation = []
     892        new_ref_pts1 = []
     893        new_ref_elev1 = []
     894        NODATA_value = -9999
     895        for i in range(6):
     896            y = (6-i)*25.0
     897            for j in range(5):
     898                x = j*25.0
     899                z = x+2*y
     900                if j == 4: z = NODATA_value # column
     901                if i == 2 and j == 2: z = NODATA_value # random
     902                if i == 5 and j == 1: z = NODATA_value
     903                if i == 1: z = NODATA_value # row
     904                if i == 3 and j == 1: z = NODATA_value # two pts/row
     905                if i == 3 and j == 3: z = NODATA_value
     906
     907                if z <> NODATA_value:
     908                    new_ref_elev1.append(z)
     909                    new_ref_pts1.append( [x,y] )
     910                   
     911                ref_points.append( [x,y] )
     912                ref_elevation.append(z)
     913                fid.write('%f ' %z)
     914            fid.write('\n')
     915
     916        fid.close()
     917
     918        #Write prj file with metadata
     919        metafilename = root+'.prj'
     920        fid = open(metafilename, 'w')
     921
     922
     923        fid.write("""Projection UTM
     924Zone 56
     925Datum WGS84
     926Zunits NO
     927Units METERS
     928Spheroid WGS84
     929Xshift 0.0000000000
     930Yshift 10000000.0000000000
     931Parameters
     932""")
     933        fid.close()
     934       
     935        #Convert to NetCDF pts
     936        convert_dem_from_ascii2netcdf(root)
     937        dem2pts(root, easting_min=2010.0, easting_max=2110.0,
     938                northing_min=3035.0, northing_max=3125.5)
     939
     940        #Check contents
     941        #Get NetCDF
     942        fid = NetCDFFile(root+'.pts', 'r')
     943
     944        # Get the variables
     945        #print fid.variables.keys()
     946        points = fid.variables['points']
     947        elevation = fid.variables['elevation']
     948
     949        #Check values
     950        assert fid.xllcorner[0] == 2010.0
     951        assert fid.yllcorner[0] == 3035.0
     952
     953        #create new reference points
     954        ref_points = []
     955        ref_elevation = []
     956        new_ref_pts2 = []
     957        new_ref_elev2 = []
     958        for i in range(4):
     959            y = (4-i)*25.0 + 25.0
     960            y_new = y + 3000.5 - 3035.0
     961            for j in range(4):
     962                x = j*25.0 + 25.0
     963                x_new = x + 2000.5 - 2010.0
     964                z = x+2*y
     965
     966                if j == 3: z = NODATA_value # column
     967                if i == 1 and j == 1: z = NODATA_value # random
     968                if i == 4 and j == 0: z = NODATA_value
     969                if i == 0: z = NODATA_value # row
     970                if i == 2 and j == 0: z = NODATA_value # two pts/row
     971                if i == 2 and j == 2: z = NODATA_value
     972
     973                if z <> NODATA_value:
     974                    new_ref_elev2.append(z)
     975                    new_ref_pts2.append( [x_new,y_new] )
     976
     977                   
     978                ref_points.append( [x_new,y_new] )
     979                ref_elevation.append(z)
     980
     981        #print points[:]
     982        #print ref_points
     983        #assert allclose(points, ref_points)
     984
     985        #print attributes[:]
     986        #print ref_elevation
     987        #assert allclose(elevation, ref_elevation)
     988
     989       
     990        assert len(points) == len(new_ref_pts2), 'length of returned points not correct'
     991        assert allclose(points, new_ref_pts2), 'points do not align'
     992
     993        #print 'elevation', elevation[:]
     994        assert len(elevation) == len(new_ref_elev2), 'length of returned elevation not correct'
     995        assert allclose(elevation, new_ref_elev2), 'elevations do not align'
     996        #Cleanup
     997        fid.close()
     998
     999
     1000        os.remove(root + '.pts')
     1001        os.remove(root + '.dem')
     1002        os.remove(root + '.asc')
     1003        os.remove(root + '.prj')
    8691004
    8701005
Note: See TracChangeset for help on using the changeset viewer.