| 1 | import os |
|---|
| 2 | import numpy as num |
|---|
| 3 | from Scientific.IO.NetCDF import NetCDFFile |
|---|
| 4 | import pylab as P |
|---|
| 5 | |
|---|
| 6 | import anuga |
|---|
| 7 | from mesh_factory import rectangular |
|---|
| 8 | from anuga.shallow_water.shallow_water_domain import Domain |
|---|
| 9 | from boundaries import Reflective_boundary |
|---|
| 10 | from anuga.coordinate_transforms.geo_reference import Geo_reference |
|---|
| 11 | from forcing import * |
|---|
| 12 | from anuga.file.sww import Write_sww |
|---|
| 13 | from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a, \ |
|---|
| 14 | netcdf_float |
|---|
| 15 | |
|---|
| 16 | def sts2sww_mesh(basename_in, basename_out=None, |
|---|
| 17 | spatial_thinning=1, verbose=False): |
|---|
| 18 | from anuga.pmesh.mesh import Mesh |
|---|
| 19 | if verbose: |
|---|
| 20 | print "Starting sts2sww_mesh" |
|---|
| 21 | |
|---|
| 22 | mean_stage=0. |
|---|
| 23 | zscale=1. |
|---|
| 24 | |
|---|
| 25 | if (basename_in[:-4]=='.sts'): |
|---|
| 26 | stsname = basename_in |
|---|
| 27 | else: |
|---|
| 28 | stsname = basename_in + '.sts' |
|---|
| 29 | |
|---|
| 30 | if verbose: print "Reading sts NetCDF file: %s" %stsname |
|---|
| 31 | infile = NetCDFFile(stsname, 'rl') |
|---|
| 32 | cellsize = infile.cellsize |
|---|
| 33 | ncols = infile.ncols |
|---|
| 34 | nrows = infile.nrows |
|---|
| 35 | no_data = infile.no_data |
|---|
| 36 | refzone = infile.zone |
|---|
| 37 | x_origin = infile.xllcorner |
|---|
| 38 | y_origin = infile.yllcorner |
|---|
| 39 | origin = num.array([x_origin[0], y_origin[0]]) |
|---|
| 40 | x = infile.variables['x'][:] |
|---|
| 41 | y = infile.variables['y'][:] |
|---|
| 42 | times = infile.variables['time'][:] |
|---|
| 43 | wind_speed_full = infile.variables['wind_speed'][:] |
|---|
| 44 | wind_angle_full = infile.variables['wind_angle'][:] |
|---|
| 45 | pressure_full = infile.variables['barometric_pressure'][:] |
|---|
| 46 | infile.close() |
|---|
| 47 | |
|---|
| 48 | number_of_points = nrows*ncols |
|---|
| 49 | points_utm = num.empty((number_of_points,2),float) |
|---|
| 50 | points_utm[:,0]=x+x_origin |
|---|
| 51 | points_utm[:,1]=y+y_origin |
|---|
| 52 | |
|---|
| 53 | thinned_indices=[] |
|---|
| 54 | for i in range(number_of_points): |
|---|
| 55 | if (i/ncols==0 or i/ncols==ncols-1 or (i/ncols)%(spatial_thinning)==0): |
|---|
| 56 | if ( i%(spatial_thinning)==0 or i%nrows==0 or i%nrows==nrows-1 ): |
|---|
| 57 | thinned_indices.append(i) |
|---|
| 58 | |
|---|
| 59 | #Spatial thinning |
|---|
| 60 | points_utm=points_utm[thinned_indices] |
|---|
| 61 | number_of_points = points_utm.shape[0] |
|---|
| 62 | number_of_timesteps = wind_speed_full.shape[0] |
|---|
| 63 | wind_speed = num.empty((number_of_timesteps,number_of_points),dtype=float) |
|---|
| 64 | wind_angle = num.empty((number_of_timesteps,number_of_points),dtype=float) |
|---|
| 65 | barometric_pressure = num.empty((number_of_timesteps,number_of_points),dtype=float) |
|---|
| 66 | if verbose: |
|---|
| 67 | print "Total number of points: ", nrows*ncols |
|---|
| 68 | print "Number of thinned points: ", number_of_points |
|---|
| 69 | for i in xrange(number_of_timesteps): |
|---|
| 70 | wind_speed[i] = wind_speed_full[i,thinned_indices] |
|---|
| 71 | wind_angle[i] = wind_angle_full[i,thinned_indices] |
|---|
| 72 | barometric_pressure[i] = pressure_full[i,thinned_indices] |
|---|
| 73 | |
|---|
| 74 | #P.plot(points_utm[:,0],points_utm[:,1],'ro') |
|---|
| 75 | #P.show() |
|---|
| 76 | |
|---|
| 77 | if verbose: |
|---|
| 78 | print "Generating sww triangulation of gems data" |
|---|
| 79 | |
|---|
| 80 | mesh = Mesh() |
|---|
| 81 | mesh.add_vertices(points_utm) |
|---|
| 82 | mesh.auto_segment(smooth_indents=True, expand_pinch=True) |
|---|
| 83 | mesh.auto_segment(mesh.shape.get_alpha() * 1.1) |
|---|
| 84 | try: |
|---|
| 85 | mesh.generate_mesh(minimum_triangle_angle=0.0, verbose=False) |
|---|
| 86 | except NoTrianglesError: |
|---|
| 87 | # This is a bit of a hack, going in and changing the data structure. |
|---|
| 88 | mesh.holes = [] |
|---|
| 89 | mesh.generate_mesh(minimum_triangle_angle=0.0, verbose=False) |
|---|
| 90 | |
|---|
| 91 | mesh_dic = mesh.Mesh2MeshList() |
|---|
| 92 | |
|---|
| 93 | points_utm=ensure_numeric(points_utm) |
|---|
| 94 | assert num.alltrue(ensure_numeric(mesh_dic['generatedpointlist']) |
|---|
| 95 | == ensure_numeric(points_utm)) |
|---|
| 96 | |
|---|
| 97 | volumes = mesh_dic['generatedtrianglelist'] |
|---|
| 98 | |
|---|
| 99 | # Write sww intro and grid stuff. |
|---|
| 100 | if (basename_out is not None and basename_out[:-4]=='.sww'): |
|---|
| 101 | swwname = basename_out |
|---|
| 102 | else: |
|---|
| 103 | swwname = basename_in + '.sww' |
|---|
| 104 | |
|---|
| 105 | if verbose: 'Output to %s' % swwname |
|---|
| 106 | |
|---|
| 107 | if verbose: |
|---|
| 108 | print "Writing sww wind and pressure field file" |
|---|
| 109 | outfile = NetCDFFile(swwname, 'wl') |
|---|
| 110 | sww = Write_sww([], ['wind_speed','wind_angle','barometric_pressure']) |
|---|
| 111 | sww.store_header(outfile, times, len(volumes), len(points_utm), |
|---|
| 112 | verbose=verbose, sww_precision='d') |
|---|
| 113 | outfile.mean_stage = mean_stage |
|---|
| 114 | outfile.zscale = zscale |
|---|
| 115 | sww.store_triangulation(outfile, points_utm, volumes, |
|---|
| 116 | refzone, |
|---|
| 117 | new_origin=origin, #check effect of this line |
|---|
| 118 | verbose=verbose) |
|---|
| 119 | |
|---|
| 120 | if verbose: |
|---|
| 121 | print 'Converting quantities' |
|---|
| 122 | |
|---|
| 123 | # Read in a time slice from the sts file and write it to the SWW file |
|---|
| 124 | |
|---|
| 125 | #print wind_angle[0,:10] |
|---|
| 126 | for i in range(len(times)): |
|---|
| 127 | sww.store_quantities(outfile, |
|---|
| 128 | slice_index=i, |
|---|
| 129 | verbose=verbose, |
|---|
| 130 | wind_speed=wind_speed[i,:], |
|---|
| 131 | wind_angle=wind_angle[i,:], |
|---|
| 132 | barometric_pressure=barometric_pressure[i,:], |
|---|
| 133 | sww_precision=num.float) |
|---|
| 134 | |
|---|
| 135 | if verbose: |
|---|
| 136 | sww.verbose_quantities(outfile) |
|---|
| 137 | outfile.close() |
|---|
