Changeset 1140 for inundation/ga


Ignore:
Timestamp:
Mar 24, 2005, 1:17:30 PM (20 years ago)
Author:
ole
Message:

Did sww2asc and tests

Location:
inundation/ga/storm_surge/pyvolution
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • inundation/ga/storm_surge/pyvolution/data_manager.py

    r1137 r1140  
    218218
    219219        self.precision = Float32 #Use single precision
    220         if hasattr(domain,'max_size'):
    221             self.max_size = domain.max_size#file size max is 2Gig
     220        if hasattr(domain, 'max_size'):
     221            self.max_size = domain.max_size #file size max is 2Gig
    222222        else:
    223223            self.max_size = max_size
    … …  
    246246            #Reference point
    247247            #Start time in seconds since the epoch (midnight 1/1/1970)
     248            #FIXME: Use Georef
    248249            fid.starttime = domain.starttime
    249250            fid.xllcorner = domain.xllcorner
    250251            fid.yllcorner = domain.yllcorner
    251             fid.zone = str(domain.zone) #FIXME: ?
     252            fid.zone = domain.zone
    252253
    253254            # dimension definitions
    … …  
    10671068   
    10681069    """
    1069     from Numeric import array, Float
     1070    from Numeric import array, Float, concatenate, NewAxis, zeros
    10701071
    10711072    #FIXME: Should be variable
    10721073    datum = 'WGS84'
    10731074    false_easting = 500000
    1074     false_northing = 10000000   
     1075    false_northing = 10000000
     1076    NODATA_value = -9999
     1077    #
     1078
    10751079   
    10761080    if quantity is None:
    … …  
    10961100    x = fid.variables['x'][:]
    10971101    y = fid.variables['y'][:]
     1102    volumes = fid.variables['volumes'][:]
     1103
     1104    vertex_points = concatenate ((x[:, NewAxis] ,y[:, NewAxis]), axis = 1)
     1105    assert len(vertex_points.shape) == 2
    10981106
    10991107    ymin = min(y); ymax = max(y)
    … …  
    11131121
    11141122    #Get quantity and reduce if applicable
    1115     q = fid.variables[quantity]
    1116 
     1123    q = fid.variables[quantity][:]
     1124
     1125   
    11171126    if len(q.shape) == 2:
    1118         q_reduced = array( number_of_points, Float )
    1119        
     1127        q_reduced = zeros( number_of_points, Float )
     1128
    11201129        for k in range(number_of_points):
    1121             q_reduced[k] = reduction( q[:,k] ) 
     1130            q_reduced[k] = reduction( q[:,k] )
    11221131
    11231132        q = q_reduced
    … …  
    11331142    prjid.write('Units         METERS\n')
    11341143    prjid.write('Spheroid      %s\n' %datum)
    1135     prjid.write('Xshift        %f\n' %false_easting)
    1136     prjid.write('Yshift        %f\n' %false_northing)
     1144    prjid.write('Xshift        %d\n' %false_easting)
     1145    prjid.write('Yshift        %d\n' %false_northing)
    11371146    prjid.write('Parameters\n')
    11381147    prjid.close()
    … …  
    11431152
    11441153    from Numeric import zeros, Float
    1145     points = zeros ( (ncols*nrows, 2), Float )
    1146 
    1147 
    1148     for i in xrange(ncols):
    1149         xg = i*cellsize
    1150         for j in xrange(nrows):
    1151             yg = j*cellsize
    1152             k = i*nrows + j
     1154    grid_points = zeros ( (ncols*nrows, 2), Float )
     1155
     1156
     1157    for i in xrange(nrows):
     1158        yg = i*cellsize       
     1159        for j in xrange(ncols):
     1160            xg = j*cellsize
     1161            k = i*ncols + j
    11531162
    11541163            #print k, xg, yg
    1155             points[k,0] = xg
    1156             points[k,1] = yg
     1164            grid_points[k,0] = xg
     1165            grid_points[k,1] = yg
    11571166
    11581167    #Interpolate
    11591168    from least_squares import Interpolation
     1169    from util import inside_polygon
    11601170   
    1161 
    1162 
     1171    interp = Interpolation(vertex_points, volumes, grid_points, alpha=0.0,
     1172                           precrop = False)
     1173
     1174    #Interpolate using quantity values
     1175    grid_values = interp.interpolate(q).flat
     1176
     1177    #Write
     1178    ascid = open(ascfile, 'w')
     1179
     1180    ascid.write('ncols         %d\n' %ncols)
     1181    ascid.write('nrows         %d\n' %nrows)
     1182    ascid.write('xllcorner     %d\n' %xllcorner)
     1183    ascid.write('yllcorner     %d\n' %yllcorner)
     1184    ascid.write('cellsize      %f\n' %cellsize)
     1185    ascid.write('NODATA_value  %d\n' %NODATA_value)                   
     1186
     1187
     1188    #Get bounding polygon from mesh
     1189    P = interp.mesh.get_boundary_polygon()
     1190
     1191    #print grid_points
     1192    #print grid_values
    11631193   
     1194    for i in range(nrows):
     1195        for j in range(ncols):
     1196            index = (nrows-i-1)*ncols+j
     1197           
     1198            if inside_polygon(grid_points[index], P):
     1199                ascid.write('%f ' %grid_values[index])               
     1200            else:
     1201                ascid.write('%d ' %NODATA_value)
     1202               
     1203        ascid.write('\n')
     1204       
    11641205    #Close
     1206    ascid.close()
    11651207    fid.close()
    11661208   

  • inundation/ga/storm_surge/pyvolution/domain.py

    r1129 r1140  
    1616    def __init__(self, coordinates, vertices, boundary = None,
    1717                 conserved_quantities = None, other_quantities = None,
    18                  tagged_elements = None, zone=None, xllcorner=0, yllcorner=0):
     18                 tagged_elements = None, zone=0, xllcorner=0, yllcorner=0):
    1919
    2020        Mesh.__init__(self, coordinates, vertices, boundary, tagged_elements)

  • inundation/ga/storm_surge/pyvolution/least_squares.py

    r1137 r1140  
    3131class ShapeError(exceptions.Exception): pass
    3232
    33 from general_mesh import General_mesh
     33#from general_mesh import General_mesh
    3434from Numeric import zeros, array, Float, Int, dot, transpose, concatenate, ArrayType
    3535from mesh import Mesh

  • inundation/ga/storm_surge/pyvolution/test_data_manager.py

    r1137 r1140  
    583583
    584584
    585     def test_sww2asc(self):
     585    def test_sww2asc_elevation(self):
    586586        """Test that sww information can be converted correctly to asc/prj
    587587        format readable by e.g. ArcView
    … …  
    593593
    594594        #Setup
    595         self.domain.filename = 'datatest' + str(id(self))
     595        self.domain.filename = 'datatest'
     596       
     597        prjfile = self.domain.filename + '.prj'
     598        ascfile = self.domain.filename + '.asc'       
     599        swwfile = self.domain.filename + '.sww'
     600       
    596601        self.domain.set_datadir('.')
    597602        self.domain.format = 'sww'
    598603        self.domain.smooth = True
    599604        self.domain.set_quantity('elevation', lambda x,y: -x-y)
     605
     606        self.domain.xllcorner = 308500
     607        self.domain.yllcorner = 6189000
     608        self.domain.zone = 56
     609       
    600610       
    601611        sww = get_dataobject(self.domain)
    … …  
    606616        sww.store_timestep('stage')
    607617
    608 
     618        cellsize = 0.25
    609619        #Check contents
    610620        #Get NetCDF
    … …  
    612622        fid = NetCDFFile(sww.filename, 'r')
    613623
    614 
    615624        # Get the variables
    616         x = fid.variables['x']
    617         y = fid.variables['y']
    618         z = fid.variables['elevation']
    619         time = fid.variables['time']
    620         stage = fid.variables['stage']
    621         #print x[:]
    622         #print y[:]
    623         #print z[:]
    624         #print stage[:]       
     625        x = fid.variables['x'][:]
     626        y = fid.variables['y'][:]
     627        z = fid.variables['elevation'][:]
     628        time = fid.variables['time'][:]
     629        stage = fid.variables['stage'][:]
     630
    625631
    626632        #Export to ascii/prj files
    627633        sww2asc(self.domain.filename,
    628634                quantity = 'elevation',                         
    629                 cellsize = 0.25,
    630                 origin = (56,0,0))
    631 
    632 
    633         raise 'Under construction (Ole)'
    634 
    635         #Check values
    636         #Q = self.domain.quantities['stage']
    637         #Q0 = Q.vertex_values[:,0]
    638         #Q1 = Q.vertex_values[:,1]
    639         #Q2 = Q.vertex_values[:,2]
    640         #
    641         #A = stage[1,:]
    642         #assert allclose(A[0], min(Q2[0], Q1[1]))
    643         #assert allclose(A[1], min(Q0[1], Q1[3], Q2[2]))
    644         #assert allclose(A[2], Q0[3])
    645         #assert allclose(A[3], min(Q0[0], Q1[5], Q2[4]))
    646         #
    647         ##Center point
    648         #assert allclose(A[4], min(Q1[0], Q2[1], Q0[2],\
    649         #                          Q0[5], Q2[6], Q1[7]))
    650 
     635                cellsize = cellsize)
     636
     637
     638        #Check prj (meta data)
     639        prjid = open(prjfile)
     640
     641        lines = prjid.readlines()
     642        prjid.close()
     643
     644        L = lines[0].strip().split()
     645        assert L[0].strip().lower() == 'projection'
     646        assert L[1].strip().lower() == 'utm'
     647
     648        L = lines[1].strip().split()
     649        assert L[0].strip().lower() == 'zone'
     650        assert L[1].strip().lower() == '56'
     651
     652        L = lines[2].strip().split()
     653        assert L[0].strip().lower() == 'datum'
     654        assert L[1].strip().lower() == 'wgs84'
     655
     656        L = lines[3].strip().split()
     657        assert L[0].strip().lower() == 'zunits'
     658        assert L[1].strip().lower() == 'no'                       
     659
     660        L = lines[4].strip().split()
     661        assert L[0].strip().lower() == 'units'
     662        assert L[1].strip().lower() == 'meters'               
     663
     664        L = lines[5].strip().split()
     665        assert L[0].strip().lower() == 'spheroid'
     666        assert L[1].strip().lower() == 'wgs84'
     667
     668        L = lines[6].strip().split()
     669        assert L[0].strip().lower() == 'xshift'
     670        assert L[1].strip().lower() == '500000'
     671
     672        L = lines[7].strip().split()
     673        assert L[0].strip().lower() == 'yshift'
     674        assert L[1].strip().lower() == '10000000'       
     675
     676        L = lines[8].strip().split()
     677        assert L[0].strip().lower() == 'parameters'
     678       
     679
     680        #Check asc file
     681        ascid = open(ascfile)
     682        lines = ascid.readlines()
     683        ascid.close()       
     684
     685        L = lines[0].strip().split()
     686        assert L[0].strip().lower() == 'ncols'
     687        assert L[1].strip().lower() == '5'
     688
     689        L = lines[1].strip().split()
     690        assert L[0].strip().lower() == 'nrows'
     691        assert L[1].strip().lower() == '5'       
     692
     693        L = lines[2].strip().split()
     694        assert L[0].strip().lower() == 'xllcorner'
     695        assert allclose(float(L[1].strip().lower()), 308500)
     696
     697        L = lines[3].strip().split()
     698        assert L[0].strip().lower() == 'yllcorner'
     699        assert allclose(float(L[1].strip().lower()), 6189000)
     700
     701        L = lines[4].strip().split()
     702        assert L[0].strip().lower() == 'cellsize'
     703        assert allclose(float(L[1].strip().lower()), cellsize)
     704
     705        L = lines[5].strip().split()
     706        assert L[0].strip() == 'NODATA_value'
     707        assert L[1].strip().lower() == '-9999'       
     708
     709        #Check grid values
     710        for j in range(5):
     711            L = lines[6+j].strip().split()           
     712            y = (4-j) * cellsize
     713            for i in range(5):
     714                assert allclose(float(L[i]), -i*cellsize - y)
     715       
    651716
    652717        fid.close()
    653718
    654719        #Cleanup
    655         os.remove(sww.filename)
     720        os.remove(prjfile)
     721        os.remove(ascfile)       
     722        os.remove(swwfile)
     723
     724
     725    def test_sww2asc_stage_reduction(self):
     726        """Test that sww information can be converted correctly to asc/prj
     727        format readable by e.g. ArcView
     728
     729        This tests the reduction of quantity stage using min
     730        """
     731
     732        import time, os
     733        from Numeric import array, zeros, allclose, Float, concatenate
     734        from Scientific.IO.NetCDF import NetCDFFile
     735
     736        #Setup
     737        self.domain.filename = 'datatest'
     738       
     739        prjfile = self.domain.filename + '.prj'
     740        ascfile = self.domain.filename + '.asc'       
     741        swwfile = self.domain.filename + '.sww'
     742       
     743        self.domain.set_datadir('.')
     744        self.domain.format = 'sww'
     745        self.domain.smooth = True
     746        self.domain.set_quantity('elevation', lambda x,y: -x-y)
     747
     748        self.domain.xllcorner = 308500
     749        self.domain.yllcorner = 6189000
     750        self.domain.zone = 56
     751       
     752       
     753        sww = get_dataobject(self.domain)
     754        sww.store_connectivity()
     755        sww.store_timestep('stage')
     756
     757        self.domain.evolve_to_end(finaltime = 0.01)
     758        sww.store_timestep('stage')
     759
     760        cellsize = 0.25
     761        #Check contents
     762        #Get NetCDF
     763
     764        fid = NetCDFFile(sww.filename, 'r')
     765
     766        # Get the variables
     767        x = fid.variables['x'][:]
     768        y = fid.variables['y'][:]
     769        z = fid.variables['elevation'][:]
     770        time = fid.variables['time'][:]
     771        stage = fid.variables['stage'][:]
     772
     773
     774        #Export to ascii/prj files
     775        sww2asc(self.domain.filename,
     776                quantity = 'stage',                         
     777                cellsize = cellsize,
     778                reduction = min)
     779
     780
     781        #Check asc file
     782        ascid = open(ascfile)
     783        lines = ascid.readlines()
     784        ascid.close()       
     785
     786        L = lines[0].strip().split()
     787        assert L[0].strip().lower() == 'ncols'
     788        assert L[1].strip().lower() == '5'
     789
     790        L = lines[1].strip().split()
     791        assert L[0].strip().lower() == 'nrows'
     792        assert L[1].strip().lower() == '5'       
     793
     794        L = lines[2].strip().split()
     795        assert L[0].strip().lower() == 'xllcorner'
     796        assert allclose(float(L[1].strip().lower()), 308500)
     797
     798        L = lines[3].strip().split()
     799        assert L[0].strip().lower() == 'yllcorner'
     800        assert allclose(float(L[1].strip().lower()), 6189000)
     801
     802        L = lines[4].strip().split()
     803        assert L[0].strip().lower() == 'cellsize'
     804        assert allclose(float(L[1].strip().lower()), cellsize)
     805
     806        L = lines[5].strip().split()
     807        assert L[0].strip() == 'NODATA_value'
     808        assert L[1].strip().lower() == '-9999'       
     809
     810       
     811        #Check grid values (where applicable)
     812        for j in range(5):
     813            if j%2 == 0:
     814                L = lines[6+j].strip().split()           
     815                jj = 4-j
     816                for i in range(5):
     817                    if i%2 == 0:
     818                        index = jj/2 + i/2*3
     819                        val0 = stage[0,index]
     820                        val1 = stage[1,index]
     821
     822                        #print i, j, index, ':', L[i], val0, val1
     823                        assert allclose(float(L[i]), min(val0, val1))
     824
     825
     826        fid.close()
     827
     828        #Cleanup
     829        os.remove(prjfile)
     830        os.remove(ascfile)       
     831        #os.remove(swwfile)
     832
     833
     834
     835
     836    def test_sww2asc_missing_points(self):
     837        """Test that sww information can be converted correctly to asc/prj
     838        format readable by e.g. ArcView
     839
     840        This test includes the writing of missing values
     841        """
     842       
     843        import time, os
     844        from Numeric import array, zeros, allclose, Float, concatenate
     845        from Scientific.IO.NetCDF import NetCDFFile
     846
     847        #Setup mesh not coinciding with rectangle.
     848        #This will cause missing values to occur in gridded data
     849
     850
     851        points = [                        [1.0, 1.0],
     852                              [0.5, 0.5], [1.0, 0.5],
     853                  [0.0, 0.0], [0.5, 0.0], [1.0, 0.0]]           
     854
     855        vertices = [ [4,1,3], [5,2,4], [1,4,2], [2,0,1]]
     856       
     857        #Create shallow water domain
     858        domain = Domain(points, vertices)
     859        domain.default_order=2
     860
     861
     862        #Set some field values
     863        domain.set_quantity('elevation', lambda x,y: -x-y)       
     864        domain.set_quantity('friction', 0.03)
     865
     866
     867        ######################
     868        # Boundary conditions
     869        B = Transmissive_boundary(domain)
     870        domain.set_boundary( {'exterior': B} )
     871
     872
     873        ######################
     874        #Initial condition - with jumps
     875
     876        bed = domain.quantities['elevation'].vertex_values
     877        stage = zeros(bed.shape, Float)
     878
     879        h = 0.3
     880        for i in range(stage.shape[0]):
     881            if i % 2 == 0:
     882                stage[i,:] = bed[i,:] + h
     883            else:
     884                stage[i,:] = bed[i,:]
     885
     886        domain.set_quantity('stage', stage)
     887        domain.distribute_to_vertices_and_edges()
     888
     889        domain.filename = 'datatest'
     890       
     891        prjfile = domain.filename + '.prj'
     892        ascfile = domain.filename + '.asc'       
     893        swwfile = domain.filename + '.sww'
     894       
     895        domain.set_datadir('.')
     896        domain.format = 'sww'
     897        domain.smooth = True
     898
     899
     900        domain.xllcorner = 308500
     901        domain.yllcorner = 6189000
     902        domain.zone = 56
     903       
     904       
     905        sww = get_dataobject(domain)
     906        sww.store_connectivity()
     907        sww.store_timestep('stage')
     908
     909        cellsize = 0.25
     910        #Check contents
     911        #Get NetCDF
     912
     913        fid = NetCDFFile(swwfile, 'r')
     914
     915        # Get the variables
     916        x = fid.variables['x'][:]
     917        y = fid.variables['y'][:]
     918        z = fid.variables['elevation'][:]
     919        time = fid.variables['time'][:]
     920
     921        #Export to ascii/prj files
     922        sww2asc(domain.filename,
     923                quantity = 'elevation',                         
     924                cellsize = cellsize)
     925
     926
     927        #Check asc file
     928        ascid = open(ascfile)
     929        lines = ascid.readlines()
     930        ascid.close()       
     931
     932        L = lines[0].strip().split()
     933        assert L[0].strip().lower() == 'ncols'
     934        assert L[1].strip().lower() == '5'
     935
     936        L = lines[1].strip().split()
     937        assert L[0].strip().lower() == 'nrows'
     938        assert L[1].strip().lower() == '5'       
     939
     940        L = lines[2].strip().split()
     941        assert L[0].strip().lower() == 'xllcorner'
     942        assert allclose(float(L[1].strip().lower()), 308500)
     943
     944        L = lines[3].strip().split()
     945        assert L[0].strip().lower() == 'yllcorner'
     946        assert allclose(float(L[1].strip().lower()), 6189000)
     947
     948        L = lines[4].strip().split()
     949        assert L[0].strip().lower() == 'cellsize'
     950        assert allclose(float(L[1].strip().lower()), cellsize)
     951
     952        L = lines[5].strip().split()
     953        assert L[0].strip() == 'NODATA_value'
     954        assert L[1].strip().lower() == '-9999'       
     955
     956
     957        #Check grid values
     958        for j in range(5):
     959            L = lines[6+j].strip().split()           
     960            y = (4-j) * cellsize
     961            for i in range(5):
     962                if i+j >= 4:
     963                    assert allclose(float(L[i]), -i*cellsize - y)
     964                else:
     965                    #Missing values
     966                    assert allclose(float(L[i]), -9999)
     967
     968
     969
     970        fid.close()
     971
     972        #Cleanup
     973        os.remove(prjfile)
     974        os.remove(ascfile)       
     975        os.remove(swwfile)
    656976
    657977
    … …  
    11751495#-------------------------------------------------------------
    11761496if __name__ == "__main__":
    1177     suite = unittest.makeSuite(Test_Data_Manager,'test')
     1497    suite = unittest.makeSuite(Test_Data_Manager,'test')   
     1498    #suite = unittest.makeSuite(Test_Data_Manager,'test_sww2asc_stage')
    11781499    runner = unittest.TextTestRunner()
    11791500    runner.run(suite)
Note: See TracChangeset for help on using the changeset viewer.