Changeset 2003


Ignore:
Timestamp:
Nov 4, 2005, 5:17:50 PM (18 years ago)
Author:
duncan
Message:

working on gippsland to sww

Location:
inundation/pyvolution
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • inundation/pyvolution/data_manager.py

    r1992 r2003  
    4848
    4949"""
    50 
    51 from Numeric import concatenate, array, Float
     50import os
     51
     52from Numeric import concatenate, array, Float, resize
    5253
    5354from coordinate_transforms.geo_reference import Geo_reference
     
    30563057    fid.close()
    30573058
    3058 def asc_csiro2sww(bath_file_in, verbose=False):
    3059     pass
    3060 
     3059def asc_csiro2sww(bath_dir, elevation_dir, sww_file, verbose=False):
     3060    """
     3061   
     3062    assume:
     3063    All files are in ersi asc format
     3064
     3065    4 types of information
     3066    bathymetry
     3067    elevation
     3068    u velocity
     3069    v velocity
     3070
     3071    Assume each type is in a seperate directory
     3072    assume one bath file with extention .000
     3073    """
     3074    from Scientific.IO.NetCDF import NetCDFFile
     3075    from Numeric import Float, Int, Int32, searchsorted, zeros, array
     3076    from Numeric import allclose, around
     3077       
     3078    from coordinate_transforms.redfearn import redfearn
     3079   
     3080    precision = Float # So if we want to change the precision its done here
     3081   
     3082    # go in to the bath dir and load the 000 file,
     3083    bath_files = os.listdir(bath_dir)
     3084    #print "bath_files",bath_files
     3085
     3086    #fixme: make more general?
     3087    bath_file = bath_files[0]
     3088    bath_dir_file =  bath_dir + os.sep + bath_file
     3089    bath_metadata,bath_grid =  _read_asc(bath_dir_file)
     3090    #print "bath_metadata",bath_metadata
     3091    #print "bath_grid",bath_grid
     3092
     3093    #Use the date.time of the bath file as a basis for
     3094    #the start time for other files
     3095    base_start = bath_file[-12:]
     3096    #print "base_start",base_start
     3097   
     3098    #go into the elevation dir and load the 000 file
     3099    elevation_dir_file = elevation_dir  + os.sep + 'el' + base_start
     3100    elevation_metadata,elevation_grid =  _read_asc(elevation_dir_file)
     3101    #print "elevation_dir_file",elevation_dir_file
     3102    #print "elevation_grid", elevation_grid
     3103   
     3104    assert bath_metadata == elevation_metadata
     3105
     3106
     3107   
     3108   
     3109    number_of_latitudes = bath_grid.shape[0]
     3110    number_of_longitudes = bath_grid.shape[1]
     3111    number_of_times = len(os.listdir(elevation_dir))
     3112    number_of_points = number_of_latitudes*number_of_longitudes
     3113    number_of_volumes = (number_of_latitudes-1)*(number_of_longitudes-1)*2
     3114
     3115    longitudes = [bath_metadata['xllcorner']+x*bath_metadata['cellsize'] for x in range(number_of_longitudes)]
     3116
     3117   
     3118    latitudes = [bath_metadata['yllcorner']+y*bath_metadata['cellsize'] for y in range(number_of_latitudes)]
     3119
     3120    # reverse order of lat, so the fist lat represents the first grid row
     3121    latitudes.reverse()
     3122   
     3123    #print "number_of_latitudes", number_of_latitudes
     3124    #print "number_of_longitudes", number_of_longitudes
     3125    #print "number_of_times", number_of_times
     3126    #print "latitudes", latitudes
     3127    #print "longitudes", longitudes
     3128   
     3129    ######### WRITE THE SWW FILE #############
     3130    # NetCDF file definition
     3131    outfile = NetCDFFile(sww_file, 'w')
     3132
     3133    #Create new file
     3134    outfile.institution = 'Geoscience Australia'
     3135    outfile.description = 'Converted from XXX'
     3136
     3137
     3138    #For sww compatibility
     3139    outfile.smoothing = 'Yes'
     3140    outfile.order = 1
     3141   
     3142
     3143    # dimension definitions
     3144    outfile.createDimension('number_of_volumes', number_of_volumes)
     3145
     3146    outfile.createDimension('number_of_vertices', 3)
     3147    outfile.createDimension('number_of_points', number_of_points)
     3148    outfile.createDimension('number_of_timesteps', number_of_times)
     3149
     3150    # variable definitions
     3151    outfile.createVariable('x', precision, ('number_of_points',))
     3152    outfile.createVariable('y', precision, ('number_of_points',))
     3153    outfile.createVariable('elevation', precision, ('number_of_points',))
     3154
     3155    #FIXME: Backwards compatibility
     3156    outfile.createVariable('z', precision, ('number_of_points',))
     3157    #################################
     3158
     3159    outfile.createVariable('volumes', Int, ('number_of_volumes',
     3160                                            'number_of_vertices'))
     3161
     3162    outfile.createVariable('time', precision,
     3163                           ('number_of_timesteps',))
     3164
     3165    outfile.createVariable('stage', precision,
     3166                           ('number_of_timesteps',
     3167                            'number_of_points'))
     3168
     3169    outfile.createVariable('xmomentum', precision,
     3170                           ('number_of_timesteps',
     3171                            'number_of_points'))
     3172
     3173    outfile.createVariable('ymomentum', precision,
     3174                           ('number_of_timesteps',
     3175                            'number_of_points'))
     3176
     3177    #Store
     3178    from coordinate_transforms.redfearn import redfearn
     3179    x = zeros(number_of_points, Float)  #Easting
     3180    y = zeros(number_of_points, Float)  #Northing
     3181
     3182    if verbose: print 'Making triangular grid'
     3183    #Get zone of 1st point.
     3184    refzone, _, _ = redfearn(latitudes[0],longitudes[0])
     3185   
     3186    vertices = {}
     3187    i = 0
     3188    for k, lat in enumerate(latitudes):       
     3189        for l, lon in enumerate(longitudes): 
     3190
     3191            vertices[l,k] = i
     3192
     3193            zone, easting, northing = redfearn(lat,lon)
     3194
     3195            msg = 'Zone boundary crossed at longitude =', lon
     3196            #assert zone == refzone, msg
     3197            #print '%7.2f %7.2f %8.2f %8.2f' %(lon, lat, easting, northing)
     3198            x[i] = easting
     3199            y[i] = northing
     3200            i += 1
     3201
     3202
     3203    #Construct 2 triangles per 'rectangular' element
     3204    volumes = []
     3205    for l in range(number_of_longitudes-1):    #X direction
     3206        for k in range(number_of_latitudes-1): #Y direction
     3207            v1 = vertices[l,k+1]
     3208            v2 = vertices[l,k]
     3209            v3 = vertices[l+1,k+1]
     3210            v4 = vertices[l+1,k]
     3211
     3212            volumes.append([v1,v2,v3]) #Upper element
     3213            volumes.append([v4,v3,v2]) #Lower element
     3214
     3215    volumes = array(volumes)
     3216
     3217    # FIXME - use coords
     3218    zone = refzone
     3219    xllcorner = min(x)
     3220    yllcorner = min(y)
     3221
     3222    outfile.xllcorner = xllcorner
     3223    outfile.yllcorner = yllcorner
     3224    outfile.zone = zone
     3225
     3226
     3227    z = resize(bath_grid,outfile.variables['z'][:].shape)
     3228    outfile.variables['x'][:] = x - xllcorner
     3229    outfile.variables['y'][:] = y - yllcorner
     3230    outfile.variables['z'][:] = z
     3231    #outfile.variables['elevation'][:] = z  #FIXME HACK
     3232    #outfile.variables['time'][:] = times   #Store time relative
     3233    #outfile.variables['time'][:] = [0]   #Store time relative
     3234    outfile.variables['volumes'][:] = volumes.astype(Int32) #On Opteron 64
     3235
     3236    # do this to create an ok sww file?
     3237    #outfile.variables['time'][:] = [0]   #Store time relative
     3238    #outfile.variables['stage'] = z
     3239    #outfile.variables['xmomentum'] = z
     3240    #outfile.variables['ymomentum'] = z
     3241   
     3242    outfile.close()
     3243   
    30613244def _read_asc(filename, verbose=False):
    30623245    """Read esri file from the following ASCII format (.asc)
     
    30713254    NODATA_value  -9999
    30723255    138.3698 137.4194 136.5062 135.5558 ..........
     3256
    30733257    """
    30743258
     
    30863270    yllcorner = float(lines.pop(0).split()[1].strip())
    30873271    cellsize = float(lines.pop(0).split()[1].strip())
    3088     NODATA_value = int(lines.pop(0).split()[1].strip())
     3272    NODATA_value = float(lines.pop(0).split()[1].strip())
    30893273   
    30903274    assert len(lines) == nrows
     
    31003284    grid = array(grid)
    31013285
    3102     return xllcorner, yllcorner, cellsize, NODATA_value, grid
    3103    
     3286    return {'xllcorner':xllcorner,
     3287            'yllcorner':yllcorner,
     3288            'cellsize':cellsize,
     3289            'NODATA_value':NODATA_value}, grid
     3290   
  • inundation/pyvolution/test_data_manager.py

    r1992 r2003  
    77from util import mean
    88import tempfile
     9import os
     10from Scientific.IO.NetCDF import NetCDFFile
    911
    1012from data_manager import *
     
    20152017
    20162018        #Cleanup
    2017         import os
    20182019        os.remove('small.sww')
    20192020
     
    26782679""")
    26792680        fid.close()
    2680         xllcorner, yllcorner, cellsize, NODATA_value, grid = \
     2681        bath_metadata, grid = \
    26812682                   data_manager._read_asc(filename, verbose=False)
    2682         self.failUnless(xllcorner  == 2000.5,  'Failed')
    2683         self.failUnless(yllcorner  == 3000.5,  'Failed')
    2684         self.failUnless(cellsize  == 25,  'Failed')
    2685         self.failUnless(NODATA_value  == -9999,  'Failed')
     2683        self.failUnless(bath_metadata['xllcorner']  == 2000.5,  'Failed')
     2684        self.failUnless(bath_metadata['yllcorner']  == 3000.5,  'Failed')
     2685        self.failUnless(bath_metadata['cellsize']  == 25,  'Failed')
     2686        self.failUnless(bath_metadata['NODATA_value']  == -9999,  'Failed')
    26862687        self.failUnless(grid[0][0]  == 97.921,  'Failed')
    26872688        self.failUnless(grid[3][6]  == 514.980,  'Failed')
     2689
     2690        os.remove(filename)
    26882691       
    2689 
     2692    def test_asc_csiro2sww(self):
     2693        import tempfile
     2694
     2695        bath_dir = tempfile.mkdtemp()
     2696        bath_dir_filename = bath_dir + os.sep +'ba19940524.000'
     2697        #bath_dir = 'bath_data_manager_test'
     2698        #print "os.getcwd( )",os.getcwd( )
     2699        elevation_dir =  tempfile.mkdtemp()
     2700        #elevation_dir = 'elev_expanded'
     2701        elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000'
     2702        elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001'
     2703       
     2704        fid = open(bath_dir_filename, 'w')
     2705        fid.write(""" ncols             3
     2706 nrows             2
     2707 xllcorner    148.00000
     2708 yllcorner    -38.00000
     2709 cellsize       0.25
     2710 nodata_value   -9999.0
     2711    90.000    60.000    30.0
     2712   -10.000   -20.000   -30.000
     2713""")
     2714        fid.close()
     2715       
     2716        fid = open(elevation_dir_filename1, 'w')
     2717        fid.write(""" ncols             3
     2718 nrows             2
     2719 xllcorner    148.00000
     2720 yllcorner    -38.00000
     2721 cellsize       0.25
     2722 nodata_value   -9999.0
     2723    90.000    60.000    30.0
     2724     0.000     0.000     0.000
     2725""")
     2726        fid.close()
     2727
     2728        fid = open(elevation_dir_filename2, 'w')
     2729        fid.write(""" ncols             3
     2730 nrows             2
     2731 xllcorner    148.00000
     2732 yllcorner    -38.00000
     2733 cellsize       0.25
     2734 nodata_value   -9999.0
     2735    90.000    60.000    30.0
     2736    10.000    10.000    10.000
     2737""")
     2738        fid.close()
     2739       
     2740        sww_file = 'test.sww'
     2741        asc_csiro2sww(bath_dir,elevation_dir, sww_file)
     2742
     2743        # check the sww file
     2744       
     2745        fid = NetCDFFile(sww_file, 'r')    #Open existing file for read
     2746        x = fid.variables['x'][:]
     2747        y = fid.variables['y'][:]
     2748        geo_ref = Geo_reference(NetCDFObject=fid)
     2749       
     2750        fid.close()
     2751   
     2752        #tidy up
     2753        os.remove(bath_dir_filename)
     2754        os.rmdir(bath_dir)
     2755
     2756        os.remove(elevation_dir_filename1)
     2757        os.remove(elevation_dir_filename2)
     2758        os.rmdir(elevation_dir)
     2759
     2760
     2761        # remove sww file
     2762        #co-ords
     2763        #Zone:   55   
     2764        #Easting:  587798.418  Northing: 5793713.242
     2765        #Latitude:   -38  0 ' 0.00000 ''  Longitude: 148  0 ' 0.00000 ''
     2766        #Grid Convergence:  0  36 ' 56.52 ''  Point Scale: 0.99969494
     2767
     2768        #Zone:   55   
     2769        #Easting:  587798.418  Northing: 5793713.242
     2770        #Latitude:   -38  0 ' 0.00000 ''  Longitude: 148  0 ' 0.00000 ''
     2771        #Grid Convergence:  0  36 ' 56.52 ''  Point Scale: 0.99969494
     2772
     2773       
    26902774#-------------------------------------------------------------
    26912775if __name__ == "__main__":
    2692     suite = unittest.makeSuite(Test_Data_Manager,'test')
     2776    suite = unittest.makeSuite(Test_Data_Manager,'test_asc_')
     2777    #suite = unittest.makeSuite(Test_Data_Manager,'test')
    26932778    #suite = unittest.makeSuite(Test_Data_Manager,'xxxtest')
    26942779    #suite = unittest.makeSuite(Test_Data_Manager,'test_sww2dem_boundingbox')   
Note: See TracChangeset for help on using the changeset viewer.