Changeset 2553


Ignore:
Timestamp:
Mar 16, 2006, 5:35:18 PM (19 years ago)
Author:
sexton
Message:

maybe this is it for dem2pts ...

Location:
inundation/pyvolution
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • inundation/pyvolution/data_manager.py

    r2541 r2553  
    12061206    i1_0 = 0
    12071207    j1_0 = 0
     1208    thisj = 0
     1209    thisi = 0
    12081210    for i in range(nrows):
    1209         y = (nrows-i)*cellsize + yllcorner
     1211        y = (nrows-i-1)*cellsize + yllcorner
    12101212        for j in range(ncols):
    12111213            x = j*cellsize + xllcorner
    12121214            if easting_min <= x <= easting_max and \
    12131215               northing_min <= y <= northing_max:
     1216                thisj = j
     1217                thisi = i
    12141218                if dem_elevation_r[i,j] == NODATA_value: nn += 1
    12151219
    12161220                if k == 0:
    12171221                    i1_0 = i
    1218                     j1_0 = j 
     1222                    j1_0 = j
    12191223                k += 1   
    12201224
    1221             index1 = j1_0
    1222             index2 = j
    1223             index3 = i1_0
    1224             index4 = i
    1225 
    1226 
    1227     index2 += 1
    1228     nrows2 = index4 - index3
    1229     ncols2 = index2 - index1 + 1
     1225    index1 = j1_0
     1226    index2 = thisj
    12301227   
    12311228    # dimension definitions
     
    12331230    ncols_in_bounding_box = int(round((easting_max-easting_min)/cellsize))
    12341231   
    1235     clippednopoints = nrows_in_bounding_box*ncols_in_bounding_box
     1232    clippednopoints = (thisi+1-i1_0)*(thisj+1-j1_0)
    12361233    nopoints = clippednopoints-nn
    12371234
     1235    clipped_dem_elev = dem_elevation_r[i1_0:thisi+1,j1_0:thisj+1]
     1236   
    12381237    if verbose and nn > 0:
    12391238        print 'There are %d values in the elevation' %totalnopoints
    1240         #print 'There are %d values in the clipped elevation' %clippednopoints
     1239        print 'There are %d values in the clipped elevation' %clippednopoints
    12411240        print 'There are %d NODATA_values in the clipped elevation' %nn
    12421241       
     
    12531252    elevation = outfile.variables['elevation']
    12541253
    1255    
     1254    lenv = index2-index1+1
    12561255    #Store data
    12571256    global_index = 0
    1258     for i in range(nrows):
     1257    #for i in range(nrows):
     1258    for i in range(i1_0,thisi+1,1):
    12591259        if verbose and i%((nrows+10)/10)==0:
    12601260            print 'Processing row %d of %d' %(i, nrows)
     
    12621262        lower_index = global_index
    12631263
    1264         v = dem_elevation_r[i,index1:index2]
     1264        v = dem_elevation_r[i,index1:index2+1]
    12651265        no_NODATA = sum(v == NODATA_value)
    12661266        if no_NODATA > 0:
    1267             newcols = len(v) - no_NODATA #ncols_in_bounding_box - no_NODATA
     1267            newcols = lenv - no_NODATA #ncols_in_bounding_box - no_NODATA
    12681268        else:
    1269             newcols = len(v) #ncols_in_bounding_box
    1270 
    1271         #print 'here', len(v), no_NODATA, newcols, v
    1272         #print 'here again', dem_elevation_r[i,:]
     1269            newcols = lenv #ncols_in_bounding_box
     1270
    12731271        telev = zeros(newcols, Float)
    12741272        tpoints = zeros((newcols, 2), Float)
     
    12761274        local_index = 0
    12771275       
    1278         y = (nrows-i)*cellsize + yllcorner
    1279         for j in range(ncols):
     1276        y = (nrows-i-1)*cellsize + yllcorner
     1277        #for j in range(ncols):
     1278        for j in range(j1_0,index2+1,1):
    12801279
    12811280            x = j*cellsize + xllcorner
  • inundation/pyvolution/test_data_manager.py

    r2526 r2553  
    571571       
    572572
    573 
    574 
    575     def test_dem2pts(self):
     573    def test_dem2pts_bounding_box_v2(self):
    576574        """Test conversion from dem in ascii format to native NetCDF xya format
    577575        """
    578576
    579577        import time, os
    580         from Numeric import array, zeros, allclose, Float, concatenate
     578        from Numeric import array, zeros, allclose, Float, concatenate, ones
    581579        from Scientific.IO.NetCDF import NetCDFFile
    582580
     
    586584        filename = root+'.asc'
    587585        fid = open(filename, 'w')
    588         fid.write("""ncols         5
    589 nrows         6
    590 xllcorner     2000.5
    591 yllcorner     3000.5
    592 cellsize      25
     586        fid.write("""ncols         10
     587nrows         10
     588xllcorner     2000
     589yllcorner     3000
     590cellsize      1
    593591NODATA_value  -9999
    594592""")
    595         #Create linear function
    596 
     593        #Create linear function       
    597594        ref_points = []
    598595        ref_elevation = []
    599         for i in range(6):
    600             y = (6-i)*25.0
    601             for j in range(5):
    602                 x = j*25.0
    603                 z = x+2*y
    604 
    605                 ref_points.append( [x,y] )
     596        x0 = 2000
     597        y = 3010
     598        yvec = range(10)
     599        xvec = range(10)
     600        z = -1
     601        for i in range(10):
     602            y = y - 1
     603            for j in range(10):
     604                x = x0 + xvec[j]
     605                z += 1
     606                ref_points.append ([x,y])
    606607                ref_elevation.append(z)
    607608                fid.write('%f ' %z)
     
    609610
    610611        fid.close()
     612           
     613        #print 'sending pts', ref_points
     614        #print 'sending elev', ref_elevation
    611615
    612616        #Write prj file with metadata
     
    629633        #Convert to NetCDF pts
    630634        convert_dem_from_ascii2netcdf(root)
    631         dem2pts(root)
     635        dem2pts(root, easting_min=2002.0, easting_max=2007.0,
     636                northing_min=3003.0, northing_max=3006.0)
    632637
    633638        #Check contents
     
    641646
    642647        #Check values
    643 
    644         #print points[:]
    645         #print ref_points
    646         assert allclose(points, ref_points)
     648        assert fid.xllcorner[0] == 2002.0
     649        assert fid.yllcorner[0] == 3003.0
     650
     651        #create new reference points
     652        newz = []       
     653        newz[0:5] = ref_elevation[32]*ones(6) + range(6)
     654        newz[6:11] = ref_elevation[42]*ones(6) + range(6)
     655        newz[12:17] = ref_elevation[52]*ones(6) + range(6)
     656        newz[18:23] = ref_elevation[62]*ones(6) + range(6)
     657        ref_elevation = []
     658        ref_elevation = newz
     659        ref_points = []
     660        new_ref_points = []
     661        x0 = 2002
     662        y = 3007
     663        yvec = range(4)
     664        xvec = range(6)
     665        for i in range(4):
     666            y = y - 1
     667            ynew = y - 3003.0
     668            for j in range(6):
     669                x = x0 + xvec[j]
     670                xnew = x - 2002.0
     671                ref_points.append([x,y])
     672                new_ref_points.append ([xnew,ynew])       
     673
     674        #print 'bounding box'
     675        #print 'new ref pts', new_ref_points[:]
     676        #print 'received pts', points[:]
     677       
     678        assert allclose(points, new_ref_points)
    647679
    648680        #print attributes[:]
    649         #print ref_elevation
     681        #print 'ref_elev', ref_elevation[:]
     682        #print 'received elev', elevation[:]
    650683        assert allclose(elevation, ref_elevation)
    651684
     
    660693
    661694
    662 
    663     def test_dem2pts_bounding_box(self):
     695    def test_dem2pts_bounding_box_removeNullvalues_v2(self):
    664696        """Test conversion from dem in ascii format to native NetCDF xya format
    665697        """
    666698
    667699        import time, os
    668         from Numeric import array, zeros, allclose, Float, concatenate
     700        from Numeric import array, zeros, allclose, Float, concatenate, ones
    669701        from Scientific.IO.NetCDF import NetCDFFile
    670702
     
    674706        filename = root+'.asc'
    675707        fid = open(filename, 'w')
    676         fid.write("""ncols         5
    677 nrows         6
    678 xllcorner     2000.5
    679 yllcorner     3000.5
    680 cellsize      25
     708        fid.write("""ncols         10
     709nrows         10
     710xllcorner     2000
     711yllcorner     3000
     712cellsize      1
    681713NODATA_value  -9999
    682714""")
    683         #Create linear function
    684 
     715        #Create linear function       
    685716        ref_points = []
    686717        ref_elevation = []
    687         for i in range(6):
    688             y = (6-i)*25.0
    689             for j in range(5):
    690                 x = j*25.0
    691                 z = x+2*y
    692 
    693                 ref_points.append( [x,y] )
    694                 ref_elevation.append(z)
    695                 fid.write('%f ' %z)
     718        x0 = 2000
     719        y = 3010
     720        yvec = range(10)
     721        xvec = range(10)
     722        z = range(100)
     723        NODATA_value = -9999
     724        count = -1
     725        for i in range(10):
     726            y = y - 1
     727            for j in range(10):
     728                x = x0 + xvec[j]
     729                ref_points.append ([x,y])
     730                count += 1
     731                if j == 4: z[count] = NODATA_value #row
     732                if j == 8: z[count] = NODATA_value #row
     733                if i == 9: z[count] = NODATA_value #column
     734                ref_elevation.append( z[count] )
     735                fid.write('%f ' %z[count])
    696736            fid.write('\n')
    697737
    698         fid.close()
     738        fid.close()             
     739       
     740        #print 'sending elev', ref_elevation
    699741
    700742        #Write prj file with metadata
     
    717759        #Convert to NetCDF pts
    718760        convert_dem_from_ascii2netcdf(root)
    719         dem2pts(root, easting_min=2010.0, easting_max=2110.0,
    720                 northing_min=3035.0, northing_max=3125.5)
     761        dem2pts(root, easting_min=2002.0, easting_max=2007.0,
     762                northing_min=3003.0, northing_max=3006.0)
    721763
    722764        #Check contents
     
    730772
    731773        #Check values
    732         assert fid.xllcorner[0] == 2010.0
    733         assert fid.yllcorner[0] == 3035.0
     774        assert fid.xllcorner[0] == 2002.0
     775        assert fid.yllcorner[0] == 3003.0
    734776
    735777        #create new reference points
     778        newz = []       
     779        newz[0:1] = ref_elevation[32]*ones(2) + range(2)
     780        newz[2:4] = ref_elevation[32]*ones(3) + range(3,6)
     781        newz[5:6] = ref_elevation[42]*ones(2) + range(2)
     782        newz[7:9] = ref_elevation[42]*ones(3) + range(3,6)
     783        newz[10:11] = ref_elevation[52]*ones(2) + range(2)
     784        newz[12:14] = ref_elevation[52]*ones(3) + range(3,6)
     785        newz[15:16] = ref_elevation[62]*ones(2) + range(2)
     786        newz[17:19] = ref_elevation[62]*ones(3) + range(3,6)
     787       
     788        ref_elevation = []
     789        ref_elevation = newz
    736790        ref_points = []
    737         ref_elevation = []
     791        new_ref_points = []
     792        x0 = 2002
     793        y = 3007
     794        yvec = range(4)
     795        xvec = range(6)
    738796        for i in range(4):
    739             y = (4-i)*25.0 + 25.0
    740             y_new = y + 3000.5 - 3035.0
    741             for j in range(4):
    742                 x = j*25.0 + 25.0
    743                 x_new = x + 2000.5 - 2010.0
    744                 z = x+2*y
    745 
    746                 ref_points.append( [x_new,y_new] )
    747                 ref_elevation.append(z)
    748 
    749         #print points[:]
    750         #print ref_points
    751         assert allclose(points, ref_points)
    752 
    753         #print attributes[:]
    754         #print ref_elevation
     797            y = y - 1
     798            ynew = y - 3003.0
     799            for j in range(6):
     800                x = x0 + xvec[j]
     801                xnew = x - 2002.0
     802                if j <> 2:
     803                    ref_points.append([x,y])
     804                    new_ref_points.append ([xnew,ynew])       
     805               
     806        assert allclose(points, new_ref_points)
    755807        assert allclose(elevation, ref_elevation)
    756808
     
    763815        os.remove(root + '.asc')
    764816        os.remove(root + '.prj')
    765 
    766     def test_dem2pts_remove_Nullvalues(self):
    767         """Test conversion from dem in ascii format to native NetCDF xya format
    768         """
    769 
    770         import time, os
    771         from Numeric import array, zeros, allclose, Float, concatenate
    772         from Scientific.IO.NetCDF import NetCDFFile
    773 
    774         #Write test asc file
    775         root = 'demtest'
    776 
    777         filename = root+'.asc'
    778         fid = open(filename, 'w')
    779         fid.write("""ncols         5
    780 nrows         6
    781 xllcorner     2000.5
    782 yllcorner     3000.5
    783 cellsize      25
    784 NODATA_value  -9999
    785 """)
    786         #Create linear function
    787         # ref_ will write all the values
    788         # new_ref_ will write the values except for NODATA_values
    789         ref_points = []
    790         ref_elevation = []
    791         new_ref_pts = []
    792         new_ref_elev = []
    793         NODATA_value = -9999
    794         for i in range(6):
    795             y = (6-i)*25.0
    796             for j in range(5):
    797                 x = j*25.0
    798                 z = x+2*y
    799                 if j == 4: z = NODATA_value # column
    800                 if i == 2 and j == 2: z = NODATA_value # random
    801                 if i == 5 and j == 1: z = NODATA_value
    802                 if i == 1: z = NODATA_value # row
    803                 if i == 3 and j == 1: z = NODATA_value # two pts/row
    804                 if i == 3 and j == 3: z = NODATA_value
    805 
    806                
    807                 if z <> NODATA_value:
    808                     new_ref_elev.append(z)
    809                     new_ref_pts.append( [x,y] )
    810                
    811                 ref_points.append( [x,y] )
    812                 ref_elevation.append(z)
    813 
    814                 fid.write('%f ' %z)
    815             fid.write('\n')
    816 
    817         fid.close()
    818 
    819        
    820         #Write prj file with metadata
    821         metafilename = root+'.prj'
    822         fid = open(metafilename, 'w')
    823 
    824 
    825         fid.write("""Projection UTM
    826 Zone 56
    827 Datum WGS84
    828 Zunits NO
    829 Units METERS
    830 Spheroid WGS84
    831 Xshift 0.0000000000
    832 Yshift 10000000.0000000000
    833 Parameters
    834 """)
    835         fid.close()
    836 
    837         #Convert to NetCDF pts
    838         convert_dem_from_ascii2netcdf(root)
    839         dem2pts(root)
    840 
    841         #Check contents
    842         #Get NetCDF
    843         fid = NetCDFFile(root+'.pts', 'r')
    844 
    845         # Get the variables
    846         #print fid.variables.keys()
    847         points = fid.variables['points']
    848         elevation = fid.variables['elevation']
    849 
    850         #Check values
    851         #print 'points', points[:]
    852         assert len(points) == len(new_ref_pts), 'length of returned points not correct'
    853         assert allclose(points, new_ref_pts), 'points do not align'
    854 
    855         #print 'elevation', elevation[:]
    856         assert len(elevation) == len(new_ref_elev), 'length of returned elevation not correct'
    857         assert allclose(elevation, new_ref_elev), 'elevations do not align'
    858 
    859         #Cleanup
    860         fid.close()
    861 
    862 
    863         os.remove(root + '.pts')
    864         os.remove(root + '.dem')
    865         os.remove(root + '.asc')
    866         os.remove(root + '.prj')
    867 
    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
    882 nrows         6
    883 xllcorner     2000.5
    884 yllcorner     3000.5
    885 cellsize      25
    886 NODATA_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
    924 Zone 56
    925 Datum WGS84
    926 Zunits NO
    927 Units METERS
    928 Spheroid WGS84
    929 Xshift 0.0000000000
    930 Yshift 10000000.0000000000
    931 Parameters
    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')
    1004 
    1005 
    1006817
    1007818
     
    40413852    runner.run(suite)
    40423853 
     3854
     3855
     3856
     3857    def test_dem2pts(self):
     3858        """Test conversion from dem in ascii format to native NetCDF xya format
     3859        """
     3860
     3861        import time, os
     3862        from Numeric import array, zeros, allclose, Float, concatenate
     3863        from Scientific.IO.NetCDF import NetCDFFile
     3864
     3865        #Write test asc file
     3866        root = 'demtest'
     3867
     3868        filename = root+'.asc'
     3869        fid = open(filename, 'w')
     3870        fid.write("""ncols         5
     3871nrows         6
     3872xllcorner     2000.5
     3873yllcorner     3000.5
     3874cellsize      25
     3875NODATA_value  -9999
     3876""")
     3877        #Create linear function
     3878
     3879        ref_points = []
     3880        ref_elevation = []
     3881        for i in range(6):
     3882            y = (6-i)*25.0
     3883            for j in range(5):
     3884                x = j*25.0
     3885                z = x+2*y
     3886
     3887                ref_points.append( [x,y] )
     3888                ref_elevation.append(z)
     3889                fid.write('%f ' %z)
     3890            fid.write('\n')
     3891
     3892        fid.close()
     3893
     3894        #Write prj file with metadata
     3895        metafilename = root+'.prj'
     3896        fid = open(metafilename, 'w')
     3897
     3898
     3899        fid.write("""Projection UTM
     3900Zone 56
     3901Datum WGS84
     3902Zunits NO
     3903Units METERS
     3904Spheroid WGS84
     3905Xshift 0.0000000000
     3906Yshift 10000000.0000000000
     3907Parameters
     3908""")
     3909        fid.close()
     3910
     3911        #Convert to NetCDF pts
     3912        convert_dem_from_ascii2netcdf(root)
     3913        dem2pts(root)
     3914
     3915        #Check contents
     3916        #Get NetCDF
     3917        fid = NetCDFFile(root+'.pts', 'r')
     3918
     3919        # Get the variables
     3920        #print fid.variables.keys()
     3921        points = fid.variables['points']
     3922        elevation = fid.variables['elevation']
     3923
     3924        #Check values
     3925
     3926        #print points[:]
     3927        #print ref_points
     3928        assert allclose(points, ref_points)
     3929
     3930        #print attributes[:]
     3931        #print ref_elevation
     3932        assert allclose(elevation, ref_elevation)
     3933
     3934        #Cleanup
     3935        fid.close()
     3936
     3937
     3938        os.remove(root + '.pts')
     3939        os.remove(root + '.dem')
     3940        os.remove(root + '.asc')
     3941        os.remove(root + '.prj')
     3942
     3943
     3944
     3945    def test_dem2pts_bounding_box(self):
     3946        """Test conversion from dem in ascii format to native NetCDF xya format
     3947        """
     3948
     3949        import time, os
     3950        from Numeric import array, zeros, allclose, Float, concatenate
     3951        from Scientific.IO.NetCDF import NetCDFFile
     3952
     3953        #Write test asc file
     3954        root = 'demtest'
     3955
     3956        filename = root+'.asc'
     3957        fid = open(filename, 'w')
     3958        fid.write("""ncols         5
     3959nrows         6
     3960xllcorner     2000.5
     3961yllcorner     3000.5
     3962cellsize      25
     3963NODATA_value  -9999
     3964""")
     3965        #Create linear function
     3966
     3967        ref_points = []
     3968        ref_elevation = []
     3969        for i in range(6):
     3970            y = (6-i)*25.0
     3971            for j in range(5):
     3972                x = j*25.0
     3973                z = x+2*y
     3974
     3975                ref_points.append( [x,y] )
     3976                ref_elevation.append(z)
     3977                fid.write('%f ' %z)
     3978            fid.write('\n')
     3979
     3980        fid.close()
     3981
     3982        #Write prj file with metadata
     3983        metafilename = root+'.prj'
     3984        fid = open(metafilename, 'w')
     3985
     3986
     3987        fid.write("""Projection UTM
     3988Zone 56
     3989Datum WGS84
     3990Zunits NO
     3991Units METERS
     3992Spheroid WGS84
     3993Xshift 0.0000000000
     3994Yshift 10000000.0000000000
     3995Parameters
     3996""")
     3997        fid.close()
     3998
     3999        #Convert to NetCDF pts
     4000        convert_dem_from_ascii2netcdf(root)
     4001        dem2pts(root, easting_min=2010.0, easting_max=2110.0,
     4002                northing_min=3035.0, northing_max=3125.5)
     4003
     4004        #Check contents
     4005        #Get NetCDF
     4006        fid = NetCDFFile(root+'.pts', 'r')
     4007
     4008        # Get the variables
     4009        #print fid.variables.keys()
     4010        points = fid.variables['points']
     4011        elevation = fid.variables['elevation']
     4012
     4013        #Check values
     4014        assert fid.xllcorner[0] == 2010.0
     4015        assert fid.yllcorner[0] == 3035.0
     4016
     4017        #create new reference points
     4018        ref_points = []
     4019        ref_elevation = []
     4020        for i in range(4):
     4021            y = (4-i)*25.0 + 25.0
     4022            y_new = y + 3000.5 - 3035.0
     4023            for j in range(4):
     4024                x = j*25.0 + 25.0
     4025                x_new = x + 2000.5 - 2010.0
     4026                z = x+2*y
     4027
     4028                ref_points.append( [x_new,y_new] )
     4029                ref_elevation.append(z)
     4030
     4031        #print points[:]
     4032        #print ref_points
     4033        assert allclose(points, ref_points)
     4034
     4035        #print attributes[:]
     4036        #print ref_elevation
     4037        assert allclose(elevation, ref_elevation)
     4038
     4039        #Cleanup
     4040        fid.close()
     4041
     4042
     4043        os.remove(root + '.pts')
     4044        os.remove(root + '.dem')
     4045        os.remove(root + '.asc')
     4046        os.remove(root + '.prj')
     4047
     4048
     4049
     4050    def test_dem2pts_remove_Nullvalues(self):
     4051        """Test conversion from dem in ascii format to native NetCDF xya format
     4052        """
     4053
     4054        import time, os
     4055        from Numeric import array, zeros, allclose, Float, concatenate
     4056        from Scientific.IO.NetCDF import NetCDFFile
     4057
     4058        #Write test asc file
     4059        root = 'demtest'
     4060
     4061        filename = root+'.asc'
     4062        fid = open(filename, 'w')
     4063        fid.write("""ncols         5
     4064nrows         6
     4065xllcorner     2000.5
     4066yllcorner     3000.5
     4067cellsize      25
     4068NODATA_value  -9999
     4069""")
     4070        #Create linear function
     4071        # ref_ will write all the values
     4072        # new_ref_ will write the values except for NODATA_values
     4073        ref_points = []
     4074        ref_elevation = []
     4075        new_ref_pts = []
     4076        new_ref_elev = []
     4077        NODATA_value = -9999
     4078        for i in range(6):
     4079            y = (6-i)*25.0
     4080            for j in range(5):
     4081                x = j*25.0
     4082                z = x+2*y
     4083                if j == 4: z = NODATA_value # column
     4084                if i == 2 and j == 2: z = NODATA_value # random
     4085                if i == 5 and j == 1: z = NODATA_value
     4086                if i == 1: z = NODATA_value # row
     4087                if i == 3 and j == 1: z = NODATA_value # two pts/row
     4088                if i == 3 and j == 3: z = NODATA_value
     4089
     4090               
     4091                if z <> NODATA_value:
     4092                    new_ref_elev.append(z)
     4093                    new_ref_pts.append( [x,y] )
     4094               
     4095                ref_points.append( [x,y] )
     4096                ref_elevation.append(z)
     4097
     4098                fid.write('%f ' %z)
     4099            fid.write('\n')
     4100
     4101        fid.close()
     4102
     4103       
     4104        #Write prj file with metadata
     4105        metafilename = root+'.prj'
     4106        fid = open(metafilename, 'w')
     4107
     4108
     4109        fid.write("""Projection UTM
     4110Zone 56
     4111Datum WGS84
     4112Zunits NO
     4113Units METERS
     4114Spheroid WGS84
     4115Xshift 0.0000000000
     4116Yshift 10000000.0000000000
     4117Parameters
     4118""")
     4119        fid.close()
     4120
     4121        #Convert to NetCDF pts
     4122        convert_dem_from_ascii2netcdf(root)
     4123        dem2pts(root)
     4124
     4125        #Check contents
     4126        #Get NetCDF
     4127        fid = NetCDFFile(root+'.pts', 'r')
     4128
     4129        # Get the variables
     4130        #print fid.variables.keys()
     4131        points = fid.variables['points']
     4132        elevation = fid.variables['elevation']
     4133
     4134        #Check values
     4135        #print 'points', points[:]
     4136        assert len(points) == len(new_ref_pts), 'length of returned points not correct'
     4137        assert allclose(points, new_ref_pts), 'points do not align'
     4138
     4139        #print 'elevation', elevation[:]
     4140        assert len(elevation) == len(new_ref_elev), 'length of returned elevation not correct'
     4141        assert allclose(elevation, new_ref_elev), 'elevations do not align'
     4142
     4143        #Cleanup
     4144        fid.close()
     4145
     4146
     4147        os.remove(root + '.pts')
     4148        os.remove(root + '.dem')
     4149        os.remove(root + '.asc')
     4150        os.remove(root + '.prj')
     4151
     4152    def test_dem2pts_bounding_box_Nullvalues(self):
     4153        """Test conversion from dem in ascii format to native NetCDF xya format
     4154        """
     4155
     4156        import time, os
     4157        from Numeric import array, zeros, allclose, Float, concatenate
     4158        from Scientific.IO.NetCDF import NetCDFFile
     4159
     4160        #Write test asc file
     4161        root = 'demtest'
     4162
     4163        filename = root+'.asc'
     4164        fid = open(filename, 'w')
     4165        fid.write("""ncols         5
     4166nrows         6
     4167xllcorner     2000.5
     4168yllcorner     3000.5
     4169cellsize      25
     4170NODATA_value  -9999
     4171""")
     4172        #Create linear function
     4173
     4174        ref_points = []
     4175        ref_elevation = []
     4176        new_ref_pts1 = []
     4177        new_ref_elev1 = []
     4178        NODATA_value = -9999
     4179        for i in range(6):
     4180            y = (6-i)*25.0
     4181            for j in range(5):
     4182                x = j*25.0
     4183                z = x+2*y
     4184                if j == 4: z = NODATA_value # column
     4185                if i == 2 and j == 2: z = NODATA_value # random
     4186                if i == 5 and j == 1: z = NODATA_value
     4187                if i == 1: z = NODATA_value # row
     4188                if i == 3 and j == 1: z = NODATA_value # two pts/row
     4189                if i == 3 and j == 3: z = NODATA_value
     4190
     4191                if z <> NODATA_value:
     4192                    new_ref_elev1.append(z)
     4193                    new_ref_pts1.append( [x,y] )
     4194                   
     4195                ref_points.append( [x,y] )
     4196                ref_elevation.append(z)
     4197                fid.write('%f ' %z)
     4198            fid.write('\n')
     4199
     4200        fid.close()
     4201
     4202        #Write prj file with metadata
     4203        metafilename = root+'.prj'
     4204        fid = open(metafilename, 'w')
     4205
     4206
     4207        fid.write("""Projection UTM
     4208Zone 56
     4209Datum WGS84
     4210Zunits NO
     4211Units METERS
     4212Spheroid WGS84
     4213Xshift 0.0000000000
     4214Yshift 10000000.0000000000
     4215Parameters
     4216""")
     4217        fid.close()
     4218       
     4219        #Convert to NetCDF pts
     4220        convert_dem_from_ascii2netcdf(root)
     4221        dem2pts(root, easting_min=2010.0, easting_max=2110.0,
     4222                northing_min=3035.0, northing_max=3125.5)
     4223
     4224        #Check contents
     4225        #Get NetCDF
     4226        fid = NetCDFFile(root+'.pts', 'r')
     4227
     4228        # Get the variables
     4229        #print fid.variables.keys()
     4230        points = fid.variables['points']
     4231        elevation = fid.variables['elevation']
     4232
     4233        #Check values
     4234        assert fid.xllcorner[0] == 2010.0
     4235        assert fid.yllcorner[0] == 3035.0
     4236
     4237        #create new reference points
     4238        ref_points = []
     4239        ref_elevation = []
     4240        new_ref_pts2 = []
     4241        new_ref_elev2 = []
     4242        for i in range(4):
     4243            y = (4-i)*25.0 + 25.0
     4244            y_new = y + 3000.5 - 3035.0
     4245            for j in range(4):
     4246                x = j*25.0 + 25.0
     4247                x_new = x + 2000.5 - 2010.0
     4248                z = x+2*y
     4249
     4250                if j == 3: z = NODATA_value # column
     4251                if i == 1 and j == 1: z = NODATA_value # random
     4252                if i == 4 and j == 0: z = NODATA_value
     4253                if i == 0: z = NODATA_value # row
     4254                if i == 2 and j == 0: z = NODATA_value # two pts/row
     4255                if i == 2 and j == 2: z = NODATA_value
     4256
     4257                if z <> NODATA_value:
     4258                    new_ref_elev2.append(z)
     4259                    new_ref_pts2.append( [x_new,y_new] )
     4260
     4261                   
     4262                ref_points.append( [x_new,y_new] )
     4263                ref_elevation.append(z)
     4264
     4265        #print points[:]
     4266        #print ref_points
     4267        #assert allclose(points, ref_points)
     4268
     4269        #print attributes[:]
     4270        #print ref_elevation
     4271        #assert allclose(elevation, ref_elevation)
     4272
     4273       
     4274        assert len(points) == len(new_ref_pts2), 'length of returned points not correct'
     4275        assert allclose(points, new_ref_pts2), 'points do not align'
     4276
     4277        #print 'elevation', elevation[:]
     4278        assert len(elevation) == len(new_ref_elev2), 'length of returned elevation not correct'
     4279        assert allclose(elevation, new_ref_elev2), 'elevations do not align'
     4280        #Cleanup
     4281        fid.close()
     4282
     4283
     4284        os.remove(root + '.pts')
     4285        os.remove(root + '.dem')
     4286        os.remove(root + '.asc')
     4287        os.remove(root + '.prj')
     4288
     4289
     4290
     4291
     4292
Note: See TracChangeset for help on using the changeset viewer.