[5265] | 1 | import os |
---|
| 2 | from struct import unpack |
---|
| 3 | from anuga.utilities.anuga_exceptions import ANUGAError |
---|
| 4 | from Numeric import Float, Int, zeros |
---|
| 5 | from Numeric import concatenate, array, Float, Int, Int32, resize, \ |
---|
| 6 | sometrue, searchsorted, zeros, allclose, around, reshape, \ |
---|
| 7 | transpose, sort, NewAxis, ArrayType, compress, take, arange, \ |
---|
| 8 | argmax, alltrue, shape, Float32, size |
---|
| 9 | |
---|
| 10 | from anuga.utilities.numerical_tools import ensure_numeric |
---|
| 11 | from anuga.config import max_float |
---|
| 12 | from anuga.coordinate_transforms.geo_reference import Geo_reference, \ |
---|
| 13 | write_NetCDF_georeference, ensure_geo_reference |
---|
| 14 | """ |
---|
| 15 | The name of the urs file names must be; |
---|
| 16 | [basename_in]_velocity-z-mux2 |
---|
| 17 | [basename_in]_velocity-e-mux2 |
---|
| 18 | [basename_in]_waveheight-n-mux2 |
---|
| 19 | """ |
---|
| 20 | def read_binary_mux2(file_in,verbose=False): |
---|
| 21 | """ |
---|
| 22 | Program for demuxing small format (mux2) files. Based upon |
---|
| 23 | demux_compact_tgdata.c written by C. Thomas Geoscience Austrlia 2007 |
---|
| 24 | Author: John Jakeman |
---|
| 25 | |
---|
| 26 | Unlike c code this program only reads in one file (source) |
---|
| 27 | """ |
---|
| 28 | eps=0.00001 |
---|
| 29 | |
---|
| 30 | ###########E###### |
---|
| 31 | # Open mux2 File # |
---|
| 32 | ################## |
---|
| 33 | |
---|
| 34 | if (os.access(file_in, os.F_OK) == 0): |
---|
| 35 | msg = 'File %s does not exist or is not accessible' %file_name |
---|
| 36 | raise IOError, msg |
---|
| 37 | if file_in[-4:]!="mux2": |
---|
| 38 | msg ="This program operates only on multiplexed files in mux2 format\nAt least one input file name does not end with mux2" |
---|
| 39 | raise IOError, msg |
---|
| 40 | |
---|
| 41 | mux_file = open(file_in, 'rb') |
---|
| 42 | |
---|
| 43 | ###################### |
---|
| 44 | # Read in the header # |
---|
| 45 | ###################### |
---|
| 46 | |
---|
| 47 | # Number of points/stations |
---|
| 48 | (num_sites,)= unpack('i',mux_file.read(4)) |
---|
| 49 | |
---|
| 50 | geolat=zeros(num_sites,Float) |
---|
| 51 | geolon=zeros(num_sites,Float) |
---|
| 52 | depth=zeros(num_sites,Float) |
---|
| 53 | dt_old=0.0 |
---|
| 54 | nt_old=0.0 |
---|
| 55 | for i in range(num_sites): |
---|
| 56 | # location and water depth in geographic coordinates |
---|
| 57 | #-180/180 |
---|
| 58 | #-90/90 |
---|
| 59 | (geolat[i],)= unpack('f',mux_file.read(4)) |
---|
| 60 | |
---|
| 61 | #dt, float - time step, seconds |
---|
| 62 | (geolon[i],) = unpack('f', mux_file.read(4)) |
---|
| 63 | |
---|
| 64 | (mcolat,) = unpack('f', mux_file.read(4)) |
---|
| 65 | (mcolon,) = unpack('f', mux_file.read(4)) |
---|
| 66 | (ig,) = unpack('i', mux_file.read(4)) |
---|
| 67 | (ilon,) = unpack('i', mux_file.read(4)) #grid point location |
---|
| 68 | (ilat,) = unpack('i', mux_file.read(4)) |
---|
| 69 | (depth[i],) = unpack('f', mux_file.read(4)) |
---|
| 70 | (centerlat,) = unpack('f', mux_file.read(4)) |
---|
| 71 | (centerlon,) = unpack('f', mux_file.read(4)) |
---|
| 72 | (offset,) = unpack('f', mux_file.read(4)) |
---|
| 73 | (az,) = unpack('f', mux_file.read(4)) |
---|
| 74 | (baz,) = unpack('f', mux_file.read(4)) |
---|
| 75 | (dt,) = unpack('f', mux_file.read(4)) |
---|
| 76 | (nt,) = unpack('i', mux_file.read(4)) |
---|
| 77 | for j in range(4): # identifier |
---|
| 78 | (id,) = unpack('f', mux_file.read(4)) |
---|
| 79 | |
---|
| 80 | msg = "Bad data in the mux file." |
---|
| 81 | if num_sites < 0: |
---|
| 82 | mux_file.close() |
---|
| 83 | raise ANUGAError, msg |
---|
| 84 | if dt < 0: |
---|
| 85 | mux_file.close() |
---|
| 86 | raise ANUGAError, msg |
---|
| 87 | if nt < 0: |
---|
| 88 | mux_file.close() |
---|
| 89 | raise ANUGAError, msg |
---|
| 90 | |
---|
| 91 | if (verbose): |
---|
| 92 | print "num_stations:",num_sites |
---|
| 93 | print "station number:",i |
---|
| 94 | print "geolat:",geolat[i] |
---|
| 95 | print "geolon:",geolon[i] |
---|
| 96 | print "mcolat:",mcolat |
---|
| 97 | print "mcolo:",mcolon |
---|
| 98 | print "ig:",ig |
---|
| 99 | print "ilon:",ilon |
---|
| 100 | print "ilat:",ilat |
---|
| 101 | print "depth:",depth[i] |
---|
| 102 | print "centerlat:",centerlat |
---|
| 103 | print "centerlon:",centerlon |
---|
| 104 | print "offset:",offset |
---|
| 105 | print "az:",az |
---|
| 106 | print "baz:",baz |
---|
| 107 | print "tstep",dt |
---|
| 108 | print "num_tsteps",nt |
---|
| 109 | print |
---|
| 110 | |
---|
| 111 | if i>0: |
---|
| 112 | msg="Stations have different time step size" |
---|
| 113 | assert (dt==dt_old),msg |
---|
| 114 | msg="Stations have different time series length" |
---|
| 115 | assert (nt==nt_old),msg |
---|
| 116 | dt_old=dt |
---|
| 117 | nt_old=nt |
---|
| 118 | |
---|
| 119 | ################ |
---|
| 120 | # Read in Data # |
---|
| 121 | ################ |
---|
| 122 | |
---|
| 123 | # Make an array to hold the start and stop steps for each station for source |
---|
| 124 | first_tstep=zeros(num_sites,Int) |
---|
| 125 | last_tstep=zeros(num_sites,Int) |
---|
| 126 | # Read the start and stop steps for each site into the array for source 0 |
---|
| 127 | for i in range(num_sites): |
---|
| 128 | (first_tstep[i],)=unpack('i', mux_file.read(4)) |
---|
| 129 | for i in range(num_sites): |
---|
| 130 | (last_tstep[i],)=unpack('i', mux_file.read(4)) |
---|
| 131 | |
---|
| 132 | # Compute the size of the data block |
---|
| 133 | numDataTotal=0 |
---|
| 134 | numData=zeros(num_sites,Int) |
---|
| 135 | last_output_step=0 |
---|
| 136 | for i in range(num_sites): |
---|
| 137 | numDataTotal+= last_tstep[i]-first_tstep[i]+1 |
---|
| 138 | numData[i]= last_tstep[i]-first_tstep[i]+1 |
---|
| 139 | last_output_step=max(last_output_step,last_tstep[i]) |
---|
| 140 | numDataTotal+=last_output_step |
---|
| 141 | |
---|
| 142 | msg="Size of Data Block is negative" |
---|
| 143 | assert numDataTotal>=0,msg |
---|
| 144 | |
---|
| 145 | muxData = zeros(numDataTotal,Float) |
---|
| 146 | for i in range(numDataTotal): |
---|
| 147 | (muxData[i],)=unpack('f', mux_file.read(4)) |
---|
| 148 | |
---|
| 149 | ################## |
---|
| 150 | # Store Mux data # |
---|
| 151 | ################## |
---|
| 152 | |
---|
| 153 | # Make arrays of starting and finishing time steps for the tide gauges |
---|
| 154 | # and fill them from the file |
---|
| 155 | data=zeros((num_sites,nt),Float) |
---|
| 156 | maximum=zeros(num_sites,Float) |
---|
| 157 | |
---|
| 158 | # Find when first station starts recording |
---|
| 159 | min_tstep = min(first_tstep) #not necessary |
---|
| 160 | # Find when all stations have stoped recording |
---|
| 161 | max_tstep = max(last_tstep) |
---|
| 162 | |
---|
| 163 | times=zeros(max_tstep,Float) |
---|
| 164 | beg=dt |
---|
| 165 | |
---|
| 166 | for i in range(num_sites): |
---|
| 167 | if (verbose): print 'reading site '+str(i) |
---|
| 168 | istart=-1 |
---|
| 169 | istop=-1 |
---|
| 170 | offset=0 |
---|
| 171 | # Update start and stop timesteps for this gauge |
---|
| 172 | if istart==-1: |
---|
| 173 | istart=first_tstep[i] |
---|
| 174 | else: |
---|
| 175 | istart=min(first_tstep[i],istart) |
---|
| 176 | if istop==-1: |
---|
| 177 | istop=last_tstep[i] |
---|
| 178 | else: |
---|
| 179 | istop=min(last_tstep[i],istop) |
---|
| 180 | for t in range(max_tstep): |
---|
| 181 | last_t=t |
---|
| 182 | offset+=1 |
---|
| 183 | # Skip records from earlier tide gauges |
---|
| 184 | for j in range(i): |
---|
| 185 | if (t+1>=first_tstep[j]) and (t+1<=last_tstep[j]): |
---|
| 186 | offset+=1 |
---|
| 187 | # Deal with the tide gauge at hand |
---|
| 188 | if (t+1>=first_tstep[i]) and (t+1<=last_tstep[i]): |
---|
| 189 | # Gauge is recording at this time] |
---|
| 190 | data[i][t]=muxData[offset] |
---|
| 191 | offset+=1 |
---|
| 192 | elif (t+1<first_tstep[i]): |
---|
| 193 | # Gauge has not yet started recording |
---|
| 194 | data[i][t]=0.0 |
---|
| 195 | else: |
---|
| 196 | # Gauge has finished recording |
---|
| 197 | data[i][t]=0.0 |
---|
| 198 | break |
---|
| 199 | # Skip records from later tide gauges |
---|
| 200 | for j in range(i+1,num_sites): |
---|
| 201 | if (t+1>=first_tstep[j]) and (t+1<=last_tstep[j]): |
---|
| 202 | offset+=1 |
---|
| 203 | if i==0: |
---|
| 204 | times[t]=beg+t*dt |
---|
| 205 | |
---|
| 206 | if (last_t<max_tstep-1): |
---|
| 207 | # The loop was exited early because the gauge had finished recording |
---|
| 208 | for t in range(last_t+1,max_tstep): |
---|
| 209 | data[i][t]=0.0 #pad until all stations have stopped recording |
---|
| 210 | |
---|
| 211 | ####################################### |
---|
| 212 | # Compute the maxima for each station # |
---|
| 213 | ####################################### |
---|
| 214 | |
---|
| 215 | for t in range(max_tstep): |
---|
| 216 | if (data[i][t] > eps): |
---|
| 217 | maximum[i]=max(data[i][t],maximum[i]) |
---|
| 218 | if i==0: |
---|
| 219 | times[t]=beg+t*dt |
---|
| 220 | |
---|
| 221 | msg = "Mux file corrupted. No positive height found at station "+str(i) |
---|
| 222 | assert maximum[i]>0, msg |
---|
| 223 | |
---|
| 224 | return times, geolat, geolon, -depth, data |
---|
| 225 | |
---|
| 226 | |
---|
| 227 | class Write_sts: |
---|
| 228 | |
---|
| 229 | sts_quantities = ['stage'] |
---|
| 230 | |
---|
| 231 | |
---|
| 232 | RANGE = '_range' |
---|
| 233 | EXTREMA = ':extrema' |
---|
| 234 | |
---|
| 235 | def __init__(self): |
---|
| 236 | pass |
---|
| 237 | |
---|
| 238 | def store_header(self, |
---|
| 239 | outfile, |
---|
| 240 | times, |
---|
| 241 | number_of_points, |
---|
| 242 | description='Converted from URS mux2 format', |
---|
| 243 | sts_precision=Float32, |
---|
| 244 | verbose=False): |
---|
| 245 | """ |
---|
| 246 | outfile - the name of the file that will be written |
---|
| 247 | times - A list of the time slice times OR a start time |
---|
| 248 | Note, if a list is given the info will be made relative. |
---|
| 249 | number_of_points - the number of urs gauges sites |
---|
| 250 | """ |
---|
| 251 | |
---|
| 252 | outfile.institution = 'Geoscience Australia' |
---|
| 253 | outfile.description = description |
---|
| 254 | |
---|
| 255 | try: |
---|
| 256 | revision_number = get_revision_number() |
---|
| 257 | except: |
---|
| 258 | revision_number = None |
---|
| 259 | # Allow None to be stored as a string |
---|
| 260 | outfile.revision_number = str(revision_number) |
---|
| 261 | |
---|
| 262 | # times - A list or array of the time slice times OR a start time |
---|
| 263 | # Start time in seconds since the epoch (midnight 1/1/1970) |
---|
| 264 | |
---|
| 265 | # This is being used to seperate one number from a list. |
---|
| 266 | # what it is actually doing is sorting lists from numeric arrays. |
---|
| 267 | if type(times) is list or type(times) is ArrayType: |
---|
| 268 | number_of_times = len(times) |
---|
| 269 | times = ensure_numeric(times) |
---|
| 270 | if number_of_times == 0: |
---|
| 271 | starttime = 0 |
---|
| 272 | else: |
---|
| 273 | starttime = times[0] |
---|
| 274 | times = times - starttime #Store relative times |
---|
| 275 | else: |
---|
| 276 | number_of_times = 0 |
---|
| 277 | starttime = times |
---|
| 278 | |
---|
| 279 | outfile.starttime = starttime |
---|
| 280 | |
---|
| 281 | # Dimension definitions |
---|
| 282 | outfile.createDimension('number_of_points', number_of_points) |
---|
| 283 | outfile.createDimension('number_of_timesteps', number_of_times) |
---|
| 284 | outfile.createDimension('numbers_in_range', 2) |
---|
| 285 | |
---|
| 286 | # Variable definitions |
---|
| 287 | outfile.createVariable('x', sts_precision, ('number_of_points',)) |
---|
| 288 | outfile.createVariable('y', sts_precision, ('number_of_points',)) |
---|
| 289 | outfile.createVariable('elevation', sts_precision,('number_of_points',)) |
---|
| 290 | |
---|
| 291 | q = 'elevation' |
---|
| 292 | outfile.createVariable(q+Write_sts.RANGE, sts_precision, |
---|
| 293 | ('numbers_in_range',)) |
---|
| 294 | |
---|
| 295 | # Initialise ranges with small and large sentinels. |
---|
| 296 | # If this was in pure Python we could have used None sensibly |
---|
| 297 | outfile.variables[q+Write_sts.RANGE][0] = max_float # Min |
---|
| 298 | outfile.variables[q+Write_sts.RANGE][1] = -max_float # Max |
---|
| 299 | |
---|
| 300 | outfile.createVariable('z', sts_precision, ('number_of_points',)) |
---|
| 301 | # Doing sts_precision instead of Float gives cast errors. |
---|
| 302 | outfile.createVariable('time', Float, ('number_of_timesteps',)) |
---|
| 303 | |
---|
| 304 | for q in Write_sts.sts_quantities: |
---|
| 305 | outfile.createVariable(q, sts_precision, |
---|
| 306 | ('number_of_timesteps', |
---|
| 307 | 'number_of_points')) |
---|
| 308 | outfile.createVariable(q+Write_sts.RANGE, sts_precision, |
---|
| 309 | ('numbers_in_range',)) |
---|
| 310 | # Initialise ranges with small and large sentinels. |
---|
| 311 | # If this was in pure Python we could have used None sensibly |
---|
| 312 | outfile.variables[q+Write_sts.RANGE][0] = max_float # Min |
---|
| 313 | outfile.variables[q+Write_sts.RANGE][1] = -max_float # Max |
---|
| 314 | |
---|
| 315 | if type(times) is list or type(times) is ArrayType: |
---|
| 316 | outfile.variables['time'][:] = times #Store time relative |
---|
| 317 | |
---|
| 318 | if verbose: |
---|
| 319 | print '------------------------------------------------' |
---|
| 320 | print 'Statistics:' |
---|
| 321 | print ' t in [%f, %f], len(t) == %d'\ |
---|
| 322 | %(min(times.flat), max(times.flat), len(times.flat)) |
---|
| 323 | |
---|
| 324 | def store_points(self, |
---|
| 325 | outfile, |
---|
| 326 | points_utm, |
---|
| 327 | elevation, zone=None, new_origin=None, |
---|
| 328 | points_georeference=None, verbose=False): |
---|
| 329 | |
---|
| 330 | """ |
---|
| 331 | points_utm - currently a list or array of the points in UTM. |
---|
| 332 | points_georeference - the georeference of the points_utm |
---|
| 333 | |
---|
| 334 | How about passing new_origin and current_origin. |
---|
| 335 | If you get both, do a convertion from the old to the new. |
---|
| 336 | |
---|
| 337 | If you only get new_origin, the points are absolute, |
---|
| 338 | convert to relative |
---|
| 339 | |
---|
| 340 | if you only get the current_origin the points are relative, store |
---|
| 341 | as relative. |
---|
| 342 | |
---|
| 343 | if you get no georefs create a new georef based on the minimums of |
---|
| 344 | points_utm. (Another option would be to default to absolute) |
---|
| 345 | |
---|
| 346 | Yes, and this is done in another part of the code. |
---|
| 347 | Probably geospatial. |
---|
| 348 | |
---|
| 349 | If you don't supply either geo_refs, then supply a zone. If not |
---|
| 350 | the default zone will be used. |
---|
| 351 | |
---|
| 352 | precondition: |
---|
| 353 | header has been called. |
---|
| 354 | """ |
---|
| 355 | |
---|
| 356 | number_of_points = len(points_utm) |
---|
| 357 | points_utm = array(points_utm) |
---|
| 358 | |
---|
| 359 | # given the two geo_refs and the points, do the stuff |
---|
| 360 | # described in the method header |
---|
| 361 | points_georeference = ensure_geo_reference(points_georeference) |
---|
| 362 | new_origin = ensure_geo_reference(new_origin) |
---|
| 363 | |
---|
| 364 | if new_origin is None and points_georeference is not None: |
---|
| 365 | points = points_utm |
---|
| 366 | geo_ref = points_georeference |
---|
| 367 | else: |
---|
| 368 | if new_origin is None: |
---|
| 369 | new_origin = Geo_reference(zone,min(points_utm[:,0]), |
---|
| 370 | min(points_utm[:,1])) |
---|
| 371 | points = new_origin.change_points_geo_ref(points_utm, |
---|
| 372 | points_georeference) |
---|
| 373 | geo_ref = new_origin |
---|
| 374 | |
---|
| 375 | # At this stage I need a georef and points |
---|
| 376 | # the points are relative to the georef |
---|
| 377 | geo_ref.write_NetCDF(outfile) |
---|
| 378 | |
---|
| 379 | x = points[:,0] |
---|
| 380 | y = points[:,1] |
---|
| 381 | z = outfile.variables['z'][:] |
---|
| 382 | |
---|
| 383 | if verbose: |
---|
| 384 | print '------------------------------------------------' |
---|
| 385 | print 'More Statistics:' |
---|
| 386 | print ' Extent (/lon):' |
---|
| 387 | print ' x in [%f, %f], len(lat) == %d'\ |
---|
| 388 | %(min(x), max(x), |
---|
| 389 | len(x)) |
---|
| 390 | print ' y in [%f, %f], len(lon) == %d'\ |
---|
| 391 | %(min(y), max(y), |
---|
| 392 | len(y)) |
---|
| 393 | print ' z in [%f, %f], len(z) == %d'\ |
---|
| 394 | %(min(elevation), max(elevation), |
---|
| 395 | len(elevation)) |
---|
| 396 | print 'geo_ref: ',geo_ref |
---|
| 397 | print '------------------------------------------------' |
---|
| 398 | |
---|
| 399 | #z = resize(bath_grid,outfile.variables['z'][:].shape) |
---|
| 400 | #print "points[:,0]", points[:,0] |
---|
| 401 | outfile.variables['x'][:] = points[:,0] #- geo_ref.get_xllcorner() |
---|
| 402 | outfile.variables['y'][:] = points[:,1] #- geo_ref.get_yllcorner() |
---|
| 403 | outfile.variables['z'][:] = elevation |
---|
| 404 | outfile.variables['elevation'][:] = elevation #FIXME HACK4 |
---|
| 405 | |
---|
| 406 | q = 'elevation' |
---|
| 407 | # This updates the _range values |
---|
| 408 | outfile.variables[q+Write_sts.RANGE][0] = min(elevation) |
---|
| 409 | outfile.variables[q+Write_sts.RANGE][1] = max(elevation) |
---|
| 410 | |
---|
| 411 | def store_quantities(self, outfile, sts_precision=Float32, |
---|
| 412 | slice_index=None, time=None, |
---|
| 413 | verbose=False, **quant): |
---|
| 414 | |
---|
| 415 | """ |
---|
| 416 | Write the quantity info. |
---|
| 417 | |
---|
| 418 | **quant is extra keyword arguments passed in. These must be |
---|
| 419 | the sts quantities, currently; stage. |
---|
| 420 | |
---|
| 421 | if the time array is already been built, use the slice_index |
---|
| 422 | to specify the index. |
---|
| 423 | |
---|
| 424 | Otherwise, use time to increase the time dimension |
---|
| 425 | |
---|
| 426 | Maybe make this general, but the viewer assumes these quantities, |
---|
| 427 | so maybe we don't want it general - unless the viewer is general |
---|
| 428 | |
---|
| 429 | precondition: |
---|
| 430 | triangulation and |
---|
| 431 | header have been called. |
---|
| 432 | """ |
---|
| 433 | if time is not None: |
---|
| 434 | file_time = outfile.variables['time'] |
---|
| 435 | slice_index = len(file_time) |
---|
| 436 | file_time[slice_index] = time |
---|
| 437 | |
---|
| 438 | # Write the conserved quantities from Domain. |
---|
| 439 | # Typically stage, xmomentum, ymomentum |
---|
| 440 | # other quantities will be ignored, silently. |
---|
| 441 | # Also write the ranges: stage_range |
---|
| 442 | for q in Write_sts.sts_quantities: |
---|
| 443 | if not quant.has_key(q): |
---|
| 444 | msg = 'STS file can not write quantity %s' %q |
---|
| 445 | raise NewQuantity, msg |
---|
| 446 | else: |
---|
| 447 | q_values = quant[q] |
---|
| 448 | outfile.variables[q][slice_index] = \ |
---|
| 449 | q_values.astype(sts_precision) |
---|
| 450 | |
---|
| 451 | # This updates the _range values |
---|
| 452 | q_range = outfile.variables[q+Write_sts.RANGE][:] |
---|
| 453 | q_values_min = min(q_values) |
---|
| 454 | if q_values_min < q_range[0]: |
---|
| 455 | outfile.variables[q+Write_sts.RANGE][0] = q_values_min |
---|
| 456 | q_values_max = max(q_values) |
---|
| 457 | if q_values_max > q_range[1]: |
---|
| 458 | outfile.variables[q+Write_sts.RANGE][1] = q_values_max |
---|
| 459 | |
---|
| 460 | def urs2sts(basename_in, basename_out = None, verbose = False, origin = None, |
---|
| 461 | mean_stage=0.0,zscale=1.0, |
---|
| 462 | minlat = None, maxlat = None, |
---|
| 463 | minlon = None, maxlon = None): |
---|
| 464 | """Convert URS mux2 format for wave propagation to sts format |
---|
| 465 | |
---|
| 466 | Specify only basename_in and read files of the form |
---|
| 467 | out_waveheight-z-mux2 |
---|
| 468 | |
---|
| 469 | Also convert latitude and longitude to UTM. All coordinates are |
---|
| 470 | assumed to be given in the GDA94 datum |
---|
| 471 | |
---|
| 472 | origin is a 3-tuple with geo referenced |
---|
| 473 | UTM coordinates (zone, easting, northing) |
---|
| 474 | """ |
---|
| 475 | import os |
---|
| 476 | from Scientific.IO.NetCDF import NetCDFFile |
---|
| 477 | from Numeric import Float, Int, Int32, searchsorted, zeros, array |
---|
| 478 | from Numeric import allclose, around |
---|
| 479 | |
---|
| 480 | precision = Float |
---|
| 481 | |
---|
| 482 | msg = 'Must use latitudes and longitudes for minlat, maxlon etc' |
---|
| 483 | |
---|
| 484 | if minlat != None: |
---|
| 485 | assert -90 < minlat < 90 , msg |
---|
| 486 | if maxlat != None: |
---|
| 487 | assert -90 < maxlat < 90 , msg |
---|
| 488 | if minlat != None: |
---|
| 489 | assert maxlat > minlat |
---|
| 490 | if minlon != None: |
---|
| 491 | assert -180 < minlon < 180 , msg |
---|
| 492 | if maxlon != None: |
---|
| 493 | assert -180 < maxlon < 180 , msg |
---|
| 494 | if minlon != None: |
---|
| 495 | assert maxlon > minlon |
---|
| 496 | |
---|
| 497 | if basename_out is None: |
---|
| 498 | stsname = basename_in + '.sts' |
---|
| 499 | else: |
---|
| 500 | stsname = basename_out + '.sts' |
---|
| 501 | |
---|
| 502 | if (verbose): print 'reading mux2 file' |
---|
| 503 | times_urs, latitudes_urs, longitudes_urs, elevations, urs_stage = read_binary_mux2(basename_in,verbose) |
---|
| 504 | |
---|
| 505 | if (minlat is not None) and (minlon is not None) and (maxlat is not None) and (maxlon is not None): |
---|
| 506 | latitudes = compress((latitudes_urs>=minlat)&(latitudes_urs<=maxlat)&(longitudes_urs>=minlon)&(longitudes_urs<=maxlon),latitudes_urs) |
---|
| 507 | longitudes = compress((latitudes_urs>=minlat)&(latitudes_urs<=maxlat)&(longitudes_urs>=minlon)&(longitudes_urs<=maxlon),longitudes_urs) |
---|
| 508 | times = compress((latitudes_urs>=minlat)&(latitudes_urs<=maxlat)&(longitudes_urs>=minlon)&(longitudes_urs<=maxlon),times_urs) |
---|
| 509 | |
---|
| 510 | |
---|
| 511 | number_of_points = latitudes.shape[0] |
---|
| 512 | number_of_times = times.shape[0] |
---|
| 513 | number_of_latitudes = latitudes.shape[0] |
---|
| 514 | number_of_longitudes = longitudes.shape[0] |
---|
| 515 | |
---|
| 516 | # NetCDF file definition |
---|
| 517 | outfile = NetCDFFile(stsname, 'w') |
---|
| 518 | |
---|
| 519 | description = 'Converted from URS mux2 files: %s'\ |
---|
| 520 | %(basename_in) |
---|
| 521 | |
---|
| 522 | # Create new file |
---|
| 523 | starttime = times[0] |
---|
| 524 | sts = Write_sts() |
---|
| 525 | sts.store_header(outfile, times,number_of_points, description=description, |
---|
| 526 | verbose=verbose,sts_precision=Float) |
---|
| 527 | |
---|
| 528 | # Store |
---|
| 529 | from anuga.coordinate_transforms.redfearn import redfearn |
---|
| 530 | x = zeros(number_of_points, Float) #Easting |
---|
| 531 | y = zeros(number_of_points, Float) #Northing |
---|
| 532 | |
---|
| 533 | # Check zone boundaries |
---|
| 534 | refzone, _, _ = redfearn(latitudes[0],longitudes[0]) |
---|
| 535 | |
---|
| 536 | for i in range(number_of_points): |
---|
| 537 | zone, easting, northing = redfearn(latitudes[i],longitudes[i]) |
---|
| 538 | x[i] = easting |
---|
| 539 | y[i] = northing |
---|
| 540 | #print zone,easting,northing |
---|
| 541 | |
---|
| 542 | if origin is None: |
---|
| 543 | origin = Geo_reference(refzone,min(x),min(y)) |
---|
| 544 | geo_ref = write_NetCDF_georeference(origin, outfile) |
---|
| 545 | |
---|
| 546 | z = elevations |
---|
| 547 | |
---|
| 548 | #print geo_ref.get_xllcorner() |
---|
| 549 | #print geo_ref.get_yllcorner() |
---|
| 550 | |
---|
| 551 | z = resize(z,outfile.variables['z'][:].shape) |
---|
| 552 | outfile.variables['x'][:] = x - geo_ref.get_xllcorner() |
---|
| 553 | outfile.variables['y'][:] = y - geo_ref.get_yllcorner() |
---|
| 554 | outfile.variables['z'][:] = z #FIXME HACK for bacwards compat. |
---|
| 555 | outfile.variables['elevation'][:] = z |
---|
| 556 | |
---|
| 557 | stage = outfile.variables['stage'] |
---|
| 558 | |
---|
| 559 | #print urs_stage.shape |
---|
| 560 | for j in range(len(times)): |
---|
| 561 | for i in range(number_of_points): |
---|
| 562 | w = zscale*urs_stage[i,j] + mean_stage |
---|
| 563 | stage[j,i] = w |
---|
| 564 | |
---|
| 565 | outfile.close() |
---|
| 566 | |
---|
| 567 | def read_sts(filename): |
---|
| 568 | from Scientific.IO.NetCDF import NetCDFFile |
---|
| 569 | from Numeric import asarray |
---|
| 570 | fid = NetCDFFile(filename+'.sts', 'r') #Open existing file for read |
---|
| 571 | time = fid.variables['time'] #Timesteps |
---|
| 572 | |
---|
| 573 | # Get the variables as Numeric arrays |
---|
| 574 | x = fid.variables['x'][:] #x-coordinates of vertices |
---|
| 575 | y = fid.variables['y'][:] #y-coordinates of vertices |
---|
| 576 | elevation = fid.variables['elevation'] #Elevation |
---|
| 577 | stage = fid.variables['stage'] #Water level |
---|
| 578 | |
---|
| 579 | starttime = fid.starttime[0] |
---|
| 580 | coordinates=transpose(asarray([x.tolist(),y.tolist()])) |
---|
| 581 | |
---|
| 582 | conserved_quantities = [] |
---|
| 583 | other_quantities = [] |
---|
| 584 | for quantity in fid.variables.keys(): |
---|
| 585 | dimensions = fid.variables[quantity].dimensions |
---|
| 586 | if 'number_of_timesteps' in dimensions: |
---|
| 587 | conserved_quantities.append(quantity) |
---|
| 588 | else: other_quantities.append(quantity) |
---|
| 589 | |
---|
| 590 | other_quantities.remove('x') |
---|
| 591 | other_quantities.remove('y') |
---|
| 592 | other_quantities.remove('z') |
---|
| 593 | try: |
---|
| 594 | other_quantities.remove('stage_range') |
---|
| 595 | other_quantities.remove('elevation_range') |
---|
| 596 | except: |
---|
| 597 | pass |
---|
| 598 | |
---|
| 599 | from pylab import plot,show |
---|
| 600 | plot(x/1000,y/1000,'o') |
---|
| 601 | show() |
---|
| 602 | |
---|
| 603 | |
---|
| 604 | #def plot_sts_gauges(): |
---|
| 605 | |
---|
| 606 | |
---|
| 607 | file_in = "big_out_waveheight-z-mux2" |
---|
| 608 | file_out="big" |
---|
| 609 | |
---|
| 610 | east=99 |
---|
| 611 | west=98 |
---|
| 612 | north=2 |
---|
| 613 | south=0 |
---|
| 614 | urs2sts(file_in,file_out,verbose=False,minlat=south,maxlat=north,minlon=west,maxlon=east) |
---|
| 615 | read_sts(file_out) |
---|