[7758] | 1 | |
---|
| 2 | ## |
---|
| 3 | # @brief Convert CSIRO ESRI file to an SWW boundary file. |
---|
| 4 | # @param bath_dir |
---|
| 5 | # @param elevation_dir |
---|
| 6 | # @param ucur_dir |
---|
| 7 | # @param vcur_dir |
---|
| 8 | # @param sww_file |
---|
| 9 | # @param minlat |
---|
| 10 | # @param maxlat |
---|
| 11 | # @param minlon |
---|
| 12 | # @param maxlon |
---|
| 13 | # @param zscale |
---|
| 14 | # @param mean_stage |
---|
| 15 | # @param fail_on_NaN |
---|
| 16 | # @param elevation_NaN_filler |
---|
| 17 | # @param bath_prefix |
---|
| 18 | # @param elevation_prefix |
---|
| 19 | # @param verbose |
---|
| 20 | # @note Also convert latitude and longitude to UTM. All coordinates are |
---|
| 21 | # assumed to be given in the GDA94 datum. |
---|
| 22 | def esri2sww(bath_dir, |
---|
| 23 | elevation_dir, |
---|
| 24 | ucur_dir, |
---|
| 25 | vcur_dir, |
---|
| 26 | sww_file, |
---|
| 27 | minlat=None, maxlat=None, |
---|
| 28 | minlon=None, maxlon=None, |
---|
| 29 | zscale=1, |
---|
| 30 | mean_stage=0, |
---|
| 31 | fail_on_NaN=True, |
---|
| 32 | elevation_NaN_filler=0, |
---|
| 33 | bath_prefix='ba', |
---|
| 34 | elevation_prefix='el', |
---|
| 35 | verbose=False): |
---|
| 36 | """ |
---|
| 37 | Produce an sww boundary file, from esri ascii data from CSIRO. |
---|
| 38 | |
---|
| 39 | Also convert latitude and longitude to UTM. All coordinates are |
---|
| 40 | assumed to be given in the GDA94 datum. |
---|
| 41 | |
---|
| 42 | assume: |
---|
| 43 | All files are in esri ascii format |
---|
| 44 | |
---|
| 45 | 4 types of information |
---|
| 46 | bathymetry |
---|
| 47 | elevation |
---|
| 48 | u velocity |
---|
| 49 | v velocity |
---|
| 50 | |
---|
| 51 | Assumptions |
---|
| 52 | The metadata of all the files is the same |
---|
| 53 | Each type is in a seperate directory |
---|
| 54 | One bath file with extention .000 |
---|
| 55 | The time period is less than 24hrs and uniform. |
---|
| 56 | """ |
---|
| 57 | |
---|
| 58 | from Scientific.IO.NetCDF import NetCDFFile |
---|
| 59 | |
---|
| 60 | from anuga.coordinate_transforms.redfearn import redfearn |
---|
| 61 | |
---|
| 62 | # So if we want to change the precision it's done here |
---|
| 63 | precision = netcdf_float |
---|
| 64 | |
---|
| 65 | # go in to the bath dir and load the only file, |
---|
| 66 | bath_files = os.listdir(bath_dir) |
---|
| 67 | bath_file = bath_files[0] |
---|
| 68 | bath_dir_file = bath_dir + os.sep + bath_file |
---|
| 69 | bath_metadata, bath_grid = _read_asc(bath_dir_file) |
---|
| 70 | |
---|
| 71 | #Use the date.time of the bath file as a basis for |
---|
| 72 | #the start time for other files |
---|
| 73 | base_start = bath_file[-12:] |
---|
| 74 | |
---|
| 75 | #go into the elevation dir and load the 000 file |
---|
| 76 | elevation_dir_file = elevation_dir + os.sep + elevation_prefix \ |
---|
| 77 | + base_start |
---|
| 78 | |
---|
| 79 | elevation_files = os.listdir(elevation_dir) |
---|
| 80 | ucur_files = os.listdir(ucur_dir) |
---|
| 81 | vcur_files = os.listdir(vcur_dir) |
---|
| 82 | elevation_files.sort() |
---|
| 83 | |
---|
| 84 | # the first elevation file should be the |
---|
| 85 | # file with the same base name as the bath data |
---|
| 86 | assert elevation_files[0] == 'el' + base_start |
---|
| 87 | |
---|
| 88 | number_of_latitudes = bath_grid.shape[0] |
---|
| 89 | number_of_longitudes = bath_grid.shape[1] |
---|
| 90 | number_of_volumes = (number_of_latitudes-1) * (number_of_longitudes-1) * 2 |
---|
| 91 | |
---|
| 92 | longitudes = [bath_metadata['xllcorner'] + x*bath_metadata['cellsize'] \ |
---|
| 93 | for x in range(number_of_longitudes)] |
---|
| 94 | latitudes = [bath_metadata['yllcorner'] + y*bath_metadata['cellsize'] \ |
---|
| 95 | for y in range(number_of_latitudes)] |
---|
| 96 | |
---|
| 97 | # reverse order of lat, so the first lat represents the first grid row |
---|
| 98 | latitudes.reverse() |
---|
| 99 | |
---|
| 100 | kmin, kmax, lmin, lmax = get_min_max_indices(latitudes[:],longitudes[:], |
---|
| 101 | minlat=minlat, maxlat=maxlat, |
---|
| 102 | minlon=minlon, maxlon=maxlon) |
---|
| 103 | |
---|
| 104 | bath_grid = bath_grid[kmin:kmax,lmin:lmax] |
---|
| 105 | latitudes = latitudes[kmin:kmax] |
---|
| 106 | longitudes = longitudes[lmin:lmax] |
---|
| 107 | number_of_latitudes = len(latitudes) |
---|
| 108 | number_of_longitudes = len(longitudes) |
---|
| 109 | number_of_times = len(os.listdir(elevation_dir)) |
---|
| 110 | number_of_points = number_of_latitudes * number_of_longitudes |
---|
| 111 | number_of_volumes = (number_of_latitudes-1) * (number_of_longitudes-1) * 2 |
---|
| 112 | |
---|
| 113 | #Work out the times |
---|
| 114 | if len(elevation_files) > 1: |
---|
| 115 | # Assume: The time period is less than 24hrs. |
---|
| 116 | time_period = (int(elevation_files[1][-3:]) \ |
---|
| 117 | - int(elevation_files[0][-3:])) * 60*60 |
---|
| 118 | times = [x*time_period for x in range(len(elevation_files))] |
---|
| 119 | else: |
---|
| 120 | times = [0.0] |
---|
| 121 | |
---|
| 122 | if verbose: |
---|
| 123 | log.critical('------------------------------------------------') |
---|
| 124 | log.critical('Statistics:') |
---|
| 125 | log.critical(' Extent (lat/lon):') |
---|
| 126 | log.critical(' lat in [%f, %f], len(lat) == %d' |
---|
| 127 | % (min(latitudes), max(latitudes), len(latitudes))) |
---|
| 128 | log.critical(' lon in [%f, %f], len(lon) == %d' |
---|
| 129 | % (min(longitudes), max(longitudes), len(longitudes))) |
---|
| 130 | log.critical(' t in [%f, %f], len(t) == %d' |
---|
| 131 | % (min(times), max(times), len(times))) |
---|
| 132 | |
---|
| 133 | ######### WRITE THE SWW FILE ############# |
---|
| 134 | |
---|
| 135 | # NetCDF file definition |
---|
| 136 | outfile = NetCDFFile(sww_file, netcdf_mode_w) |
---|
| 137 | |
---|
| 138 | #Create new file |
---|
| 139 | outfile.institution = 'Geoscience Australia' |
---|
| 140 | outfile.description = 'Converted from XXX' |
---|
| 141 | |
---|
| 142 | #For sww compatibility |
---|
| 143 | outfile.smoothing = 'Yes' |
---|
| 144 | outfile.order = 1 |
---|
| 145 | |
---|
| 146 | #Start time in seconds since the epoch (midnight 1/1/1970) |
---|
| 147 | outfile.starttime = starttime = times[0] |
---|
| 148 | |
---|
| 149 | # dimension definitions |
---|
| 150 | outfile.createDimension('number_of_volumes', number_of_volumes) |
---|
| 151 | outfile.createDimension('number_of_vertices', 3) |
---|
| 152 | outfile.createDimension('number_of_points', number_of_points) |
---|
| 153 | outfile.createDimension('number_of_timesteps', number_of_times) |
---|
| 154 | |
---|
| 155 | # variable definitions |
---|
| 156 | outfile.createVariable('x', precision, ('number_of_points',)) |
---|
| 157 | outfile.createVariable('y', precision, ('number_of_points',)) |
---|
| 158 | outfile.createVariable('elevation', precision, ('number_of_points',)) |
---|
| 159 | |
---|
| 160 | #FIXME: Backwards compatibility |
---|
| 161 | #outfile.createVariable('z', precision, ('number_of_points',)) |
---|
| 162 | ################################# |
---|
| 163 | |
---|
| 164 | outfile.createVariable('volumes', netcdf_int, ('number_of_volumes', |
---|
| 165 | 'number_of_vertices')) |
---|
| 166 | |
---|
| 167 | outfile.createVariable('time', precision, ('number_of_timesteps',)) |
---|
| 168 | |
---|
| 169 | outfile.createVariable('stage', precision, ('number_of_timesteps', |
---|
| 170 | 'number_of_points')) |
---|
| 171 | |
---|
| 172 | outfile.createVariable('xmomentum', precision, ('number_of_timesteps', |
---|
| 173 | 'number_of_points')) |
---|
| 174 | |
---|
| 175 | outfile.createVariable('ymomentum', precision, ('number_of_timesteps', |
---|
| 176 | 'number_of_points')) |
---|
| 177 | |
---|
| 178 | #Store |
---|
| 179 | from anuga.coordinate_transforms.redfearn import redfearn |
---|
| 180 | |
---|
| 181 | x = num.zeros(number_of_points, num.float) #Easting |
---|
| 182 | y = num.zeros(number_of_points, num.float) #Northing |
---|
| 183 | |
---|
| 184 | if verbose: log.critical('Making triangular grid') |
---|
| 185 | |
---|
| 186 | #Get zone of 1st point. |
---|
| 187 | refzone, _, _ = redfearn(latitudes[0], longitudes[0]) |
---|
| 188 | |
---|
| 189 | vertices = {} |
---|
| 190 | i = 0 |
---|
| 191 | for k, lat in enumerate(latitudes): |
---|
| 192 | for l, lon in enumerate(longitudes): |
---|
| 193 | vertices[l,k] = i |
---|
| 194 | |
---|
| 195 | zone, easting, northing = redfearn(lat, lon) |
---|
| 196 | |
---|
| 197 | #msg = 'Zone boundary crossed at longitude =', lon |
---|
| 198 | #assert zone == refzone, msg |
---|
| 199 | #print '%7.2f %7.2f %8.2f %8.2f' %(lon, lat, easting, northing) |
---|
| 200 | x[i] = easting |
---|
| 201 | y[i] = northing |
---|
| 202 | i += 1 |
---|
| 203 | |
---|
| 204 | #Construct 2 triangles per 'rectangular' element |
---|
| 205 | volumes = [] |
---|
| 206 | for l in range(number_of_longitudes-1): #X direction |
---|
| 207 | for k in range(number_of_latitudes-1): #Y direction |
---|
| 208 | v1 = vertices[l,k+1] |
---|
| 209 | v2 = vertices[l,k] |
---|
| 210 | v3 = vertices[l+1,k+1] |
---|
| 211 | v4 = vertices[l+1,k] |
---|
| 212 | |
---|
| 213 | #Note, this is different to the ferrit2sww code |
---|
| 214 | #since the order of the lats is reversed. |
---|
| 215 | volumes.append([v1,v3,v2]) #Upper element |
---|
| 216 | volumes.append([v4,v2,v3]) #Lower element |
---|
| 217 | |
---|
| 218 | volumes = num.array(volumes, num.int) #array default# |
---|
| 219 | |
---|
| 220 | geo_ref = Geo_reference(refzone, min(x), min(y)) |
---|
| 221 | geo_ref.write_NetCDF(outfile) |
---|
| 222 | |
---|
| 223 | # This will put the geo ref in the middle |
---|
| 224 | #geo_ref = Geo_reference(refzone, (max(x)+min(x))/2., (max(x)+min(y))/2.) |
---|
| 225 | |
---|
| 226 | if verbose: |
---|
| 227 | log.critical('------------------------------------------------') |
---|
| 228 | log.critical('More Statistics:') |
---|
| 229 | log.critical(' Extent (/lon):') |
---|
| 230 | log.critical(' x in [%f, %f], len(lat) == %d' |
---|
| 231 | % (min(x), max(x), len(x))) |
---|
| 232 | log.critical(' y in [%f, %f], len(lon) == %d' |
---|
| 233 | % (min(y), max(y), len(y))) |
---|
| 234 | log.critical('geo_ref: ', geo_ref) |
---|
| 235 | |
---|
| 236 | z = num.resize(bath_grid,outfile.variables['elevation'][:].shape) |
---|
| 237 | outfile.variables['x'][:] = x - geo_ref.get_xllcorner() |
---|
| 238 | outfile.variables['y'][:] = y - geo_ref.get_yllcorner() |
---|
| 239 | # FIXME (Ole): Remove once viewer has been recompiled and changed |
---|
| 240 | # to use elevation instead of z |
---|
| 241 | #outfile.variables['z'][:] = z |
---|
| 242 | outfile.variables['elevation'][:] = z |
---|
| 243 | outfile.variables['volumes'][:] = volumes.astype(num.int32) # On Opteron 64 |
---|
| 244 | |
---|
| 245 | stage = outfile.variables['stage'] |
---|
| 246 | xmomentum = outfile.variables['xmomentum'] |
---|
| 247 | ymomentum = outfile.variables['ymomentum'] |
---|
| 248 | |
---|
| 249 | outfile.variables['time'][:] = times #Store time relative |
---|
| 250 | |
---|
| 251 | if verbose: log.critical('Converting quantities') |
---|
| 252 | |
---|
| 253 | n = number_of_times |
---|
| 254 | for j in range(number_of_times): |
---|
| 255 | # load in files |
---|
| 256 | elevation_meta, elevation_grid = \ |
---|
| 257 | _read_asc(elevation_dir + os.sep + elevation_files[j]) |
---|
| 258 | |
---|
| 259 | _, u_momentum_grid = _read_asc(ucur_dir + os.sep + ucur_files[j]) |
---|
| 260 | _, v_momentum_grid = _read_asc(vcur_dir + os.sep + vcur_files[j]) |
---|
| 261 | |
---|
| 262 | #cut matrix to desired size |
---|
| 263 | elevation_grid = elevation_grid[kmin:kmax,lmin:lmax] |
---|
| 264 | u_momentum_grid = u_momentum_grid[kmin:kmax,lmin:lmax] |
---|
| 265 | v_momentum_grid = v_momentum_grid[kmin:kmax,lmin:lmax] |
---|
| 266 | |
---|
| 267 | # handle missing values |
---|
| 268 | missing = (elevation_grid == elevation_meta['NODATA_value']) |
---|
| 269 | if num.sometrue (missing): |
---|
| 270 | if fail_on_NaN: |
---|
| 271 | msg = 'File %s contains missing values' \ |
---|
| 272 | % (elevation_files[j]) |
---|
| 273 | raise DataMissingValuesError, msg |
---|
| 274 | else: |
---|
| 275 | elevation_grid = elevation_grid*(missing==0) \ |
---|
| 276 | + missing*elevation_NaN_filler |
---|
| 277 | |
---|
| 278 | if verbose and j % ((n+10)/10) == 0: log.critical(' Doing %d of %d' |
---|
| 279 | % (j, n)) |
---|
| 280 | |
---|
| 281 | i = 0 |
---|
| 282 | for k in range(number_of_latitudes): #Y direction |
---|
| 283 | for l in range(number_of_longitudes): #X direction |
---|
| 284 | w = zscale*elevation_grid[k,l] + mean_stage |
---|
| 285 | stage[j,i] = w |
---|
| 286 | h = w - z[i] |
---|
| 287 | xmomentum[j,i] = u_momentum_grid[k,l]*h |
---|
| 288 | ymomentum[j,i] = v_momentum_grid[k,l]*h |
---|
| 289 | i += 1 |
---|
| 290 | |
---|
| 291 | outfile.close() |
---|
| 292 | |
---|
| 293 | |
---|
| 294 | |
---|
| 295 | def _read_asc(filename, verbose=False): |
---|
| 296 | """Read esri file from the following ASCII format (.asc) |
---|
| 297 | |
---|
| 298 | Example: |
---|
| 299 | |
---|
| 300 | ncols 3121 |
---|
| 301 | nrows 1800 |
---|
| 302 | xllcorner 722000 |
---|
| 303 | yllcorner 5893000 |
---|
| 304 | cellsize 25 |
---|
| 305 | NODATA_value -9999 |
---|
| 306 | 138.3698 137.4194 136.5062 135.5558 .......... |
---|
| 307 | |
---|
| 308 | filename Path to the file to read. |
---|
| 309 | verbose True if this function is to be verbose. |
---|
| 310 | """ |
---|
| 311 | |
---|
| 312 | datafile = open(filename) |
---|
| 313 | |
---|
| 314 | if verbose: log.critical('Reading DEM from %s' % filename) |
---|
| 315 | |
---|
| 316 | lines = datafile.readlines() |
---|
| 317 | datafile.close() |
---|
| 318 | |
---|
| 319 | if verbose: log.critical('Got %d lines' % len(lines)) |
---|
| 320 | |
---|
| 321 | ncols = int(lines.pop(0).split()[1].strip()) |
---|
| 322 | nrows = int(lines.pop(0).split()[1].strip()) |
---|
| 323 | xllcorner = float(lines.pop(0).split()[1].strip()) |
---|
| 324 | yllcorner = float(lines.pop(0).split()[1].strip()) |
---|
| 325 | cellsize = float(lines.pop(0).split()[1].strip()) |
---|
| 326 | NODATA_value = float(lines.pop(0).split()[1].strip()) |
---|
| 327 | |
---|
| 328 | assert len(lines) == nrows |
---|
| 329 | |
---|
| 330 | #Store data |
---|
| 331 | grid = [] |
---|
| 332 | |
---|
| 333 | n = len(lines) |
---|
| 334 | for i, line in enumerate(lines): |
---|
| 335 | cells = line.split() |
---|
| 336 | assert len(cells) == ncols |
---|
| 337 | grid.append(num.array([float(x) for x in cells])) |
---|
| 338 | grid = num.array(grid) |
---|
| 339 | |
---|
| 340 | return {'xllcorner':xllcorner, |
---|
| 341 | 'yllcorner':yllcorner, |
---|
| 342 | 'cellsize':cellsize, |
---|
| 343 | 'NODATA_value':NODATA_value}, grid |
---|