Changeset 2555


Ignore:
Timestamp:
Mar 17, 2006, 2:03:19 PM (18 years ago)
Author:
ole
Message:

More testing of dem2pts (It still works :-)

File:
1 edited

Legend:

Unmodified
Added
Removed
  • inundation/pyvolution/test_data_manager.py

    r2553 r2555  
    651651        #create new reference points
    652652        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)
     653        newz[0:5] = ref_elevation[32:38]
     654        newz[6:11] = ref_elevation[42:48]
     655        newz[12:17] = ref_elevation[52:58]
     656        newz[18:23] = ref_elevation[62:68]
    657657        ref_elevation = []
    658658        ref_elevation = newz
    659659        ref_points = []
    660         new_ref_points = []
    661660        x0 = 2002
    662661        y = 3007
     
    669668                x = x0 + xvec[j]
    670669                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)
    679 
    680         #print attributes[:]
    681         #print 'ref_elev', ref_elevation[:]
    682         #print 'received elev', elevation[:]
     670                ref_points.append ([xnew,ynew]) #Relative point values       
     671
     672        assert allclose(points, ref_points)
     673
    683674        assert allclose(elevation, ref_elevation)
    684675
     
    720711        yvec = range(10)
    721712        xvec = range(10)
    722         z = range(100)
     713        #z = range(100)
     714        z = zeros(100)
    723715        NODATA_value = -9999
    724716        count = -1
     
    729721                ref_points.append ([x,y])
    730722                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
     723                z[count] = (4*i - 3*j)%13
     724                if j == 4: z[count] = NODATA_value #column inside clipping region
     725                if j == 8: z[count] = NODATA_value #column outside clipping region
     726                if i == 9: z[count] = NODATA_value #row outside clipping region
     727                if i == 4 and j == 6: z[count] = NODATA_value #arbitrary point inside clipping region               
    734728                ref_elevation.append( z[count] )
    735729                fid.write('%f ' %z[count])
     
    776770
    777771        #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 = []
     772        newz = zeros(19)
     773        newz[0:2] = ref_elevation[32:34]
     774        newz[2:5] = ref_elevation[35:38]
     775        newz[5:7] = ref_elevation[42:44]
     776        newz[7] = ref_elevation[45]
     777        newz[8] = ref_elevation[47]
     778        newz[9:11] = ref_elevation[52:54]
     779        newz[11:14] = ref_elevation[55:58]
     780        newz[14:16] = ref_elevation[62:64]
     781        newz[16:19] = ref_elevation[65:68]
     782
     783       
    789784        ref_elevation = newz
    790785        ref_points = []
     
    800795                x = x0 + xvec[j]
    801796                xnew = x - 2002.0
    802                 if j <> 2:
     797                if j <> 2 and (i<>1 or j<>4):   
    803798                    ref_points.append([x,y])
    804                     new_ref_points.append ([xnew,ynew])       
     799                    new_ref_points.append ([xnew,ynew])
     800
    805801               
    806802        assert allclose(points, new_ref_points)
    807803        assert allclose(elevation, ref_elevation)
     804
     805        #Cleanup
     806        fid.close()
     807
     808
     809        os.remove(root + '.pts')
     810        os.remove(root + '.dem')
     811        os.remove(root + '.asc')
     812        os.remove(root + '.prj')
     813
     814
     815    def test_dem2pts_bounding_box_removeNullvalues_v3(self):
     816        """Test conversion from dem in ascii format to native NetCDF xya format
     817        Check missing values on clipping boundary
     818        """
     819
     820        import time, os
     821        from Numeric import array, zeros, allclose, Float, concatenate, ones
     822        from Scientific.IO.NetCDF import NetCDFFile
     823
     824        #Write test asc file
     825        root = 'demtest'
     826
     827        filename = root+'.asc'
     828        fid = open(filename, 'w')
     829        fid.write("""ncols         10
     830nrows         10
     831xllcorner     2000
     832yllcorner     3000
     833cellsize      1
     834NODATA_value  -9999
     835""")
     836        #Create linear function       
     837        ref_points = []
     838        ref_elevation = []
     839        x0 = 2000
     840        y = 3010
     841        yvec = range(10)
     842        xvec = range(10)
     843        #z = range(100)
     844        z = zeros(100)
     845        NODATA_value = -9999
     846        count = -1
     847        for i in range(10):
     848            y = y - 1
     849            for j in range(10):
     850                x = x0 + xvec[j]
     851                ref_points.append ([x,y])
     852                count += 1
     853                z[count] = (4*i - 3*j)%13
     854                if j == 4: z[count] = NODATA_value #column inside clipping region
     855                if j == 8: z[count] = NODATA_value #column outside clipping region
     856                if i == 6: z[count] = NODATA_value #row on clipping boundary
     857                if i == 4 and j == 6: z[count] = NODATA_value #arbitrary point inside clipping region               
     858                ref_elevation.append( z[count] )
     859                fid.write('%f ' %z[count])
     860            fid.write('\n')
     861
     862        fid.close()             
     863       
     864        #print 'sending elev', ref_elevation
     865
     866        #Write prj file with metadata
     867        metafilename = root+'.prj'
     868        fid = open(metafilename, 'w')
     869
     870
     871        fid.write("""Projection UTM
     872Zone 56
     873Datum WGS84
     874Zunits NO
     875Units METERS
     876Spheroid WGS84
     877Xshift 0.0000000000
     878Yshift 10000000.0000000000
     879Parameters
     880""")
     881        fid.close()
     882
     883        #Convert to NetCDF pts
     884        convert_dem_from_ascii2netcdf(root)
     885        dem2pts(root, easting_min=2002.0, easting_max=2007.0,
     886                northing_min=3003.0, northing_max=3006.0)
     887
     888        #Check contents
     889        #Get NetCDF
     890        fid = NetCDFFile(root+'.pts', 'r')
     891
     892        # Get the variables
     893        #print fid.variables.keys()
     894        points = fid.variables['points']
     895        elevation = fid.variables['elevation']
     896
     897        #Check values
     898        assert fid.xllcorner[0] == 2002.0
     899        assert fid.yllcorner[0] == 3003.0
     900
     901        #create new reference points
     902        newz = zeros(14)
     903        newz[0:2] = ref_elevation[32:34]
     904        newz[2:5] = ref_elevation[35:38]
     905        newz[5:7] = ref_elevation[42:44]
     906        newz[7] = ref_elevation[45]
     907        newz[8] = ref_elevation[47]
     908        newz[9:11] = ref_elevation[52:54]
     909        newz[11:14] = ref_elevation[55:58]
     910
     911
     912       
     913        ref_elevation = newz
     914        ref_points = []
     915        new_ref_points = []
     916        x0 = 2002
     917        y = 3007
     918        yvec = range(4)
     919        xvec = range(6)
     920        for i in range(4):
     921            y = y - 1
     922            ynew = y - 3003.0
     923            for j in range(6):
     924                x = x0 + xvec[j]
     925                xnew = x - 2002.0
     926                if j <> 2 and (i<>1 or j<>4) and i<>3:   
     927                    ref_points.append([x,y])
     928                    new_ref_points.append ([xnew,ynew])
     929
     930
     931        #print points[:],points[:].shape
     932        #print new_ref_points, len(new_ref_points)
     933
     934        assert allclose(elevation, ref_elevation)                   
     935        assert allclose(points, new_ref_points)
     936
    808937
    809938        #Cleanup
Note: See TracChangeset for help on using the changeset viewer.