Changeset 2003
- Timestamp:
- Nov 4, 2005, 5:17:50 PM (18 years ago)
- Location:
- inundation/pyvolution
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
inundation/pyvolution/data_manager.py
r1992 r2003 48 48 49 49 """ 50 51 from Numeric import concatenate, array, Float 50 import os 51 52 from Numeric import concatenate, array, Float, resize 52 53 53 54 from coordinate_transforms.geo_reference import Geo_reference … … 3056 3057 fid.close() 3057 3058 3058 def asc_csiro2sww(bath_file_in, verbose=False): 3059 pass 3060 3059 def 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 3061 3244 def _read_asc(filename, verbose=False): 3062 3245 """Read esri file from the following ASCII format (.asc) … … 3071 3254 NODATA_value -9999 3072 3255 138.3698 137.4194 136.5062 135.5558 .......... 3256 3073 3257 """ 3074 3258 … … 3086 3270 yllcorner = float(lines.pop(0).split()[1].strip()) 3087 3271 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()) 3089 3273 3090 3274 assert len(lines) == nrows … … 3100 3284 grid = array(grid) 3101 3285 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 7 7 from util import mean 8 8 import tempfile 9 import os 10 from Scientific.IO.NetCDF import NetCDFFile 9 11 10 12 from data_manager import * … … 2015 2017 2016 2018 #Cleanup 2017 import os2018 2019 os.remove('small.sww') 2019 2020 … … 2678 2679 """) 2679 2680 fid.close() 2680 xllcorner, yllcorner, cellsize, NODATA_value, grid = \2681 bath_metadata, grid = \ 2681 2682 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') 2686 2687 self.failUnless(grid[0][0] == 97.921, 'Failed') 2687 2688 self.failUnless(grid[3][6] == 514.980, 'Failed') 2689 2690 os.remove(filename) 2688 2691 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 2690 2774 #------------------------------------------------------------- 2691 2775 if __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') 2693 2778 #suite = unittest.makeSuite(Test_Data_Manager,'xxxtest') 2694 2779 #suite = unittest.makeSuite(Test_Data_Manager,'test_sww2dem_boundingbox')
Note: See TracChangeset
for help on using the changeset viewer.