source: inundation/pyvolution/data_manager_obsolete_stuff.py @ 3290

Last change on this file since 3290 was 2344, checked in by ole, 19 years ago

Moved obsolete stuff out from data_manager
File size: 10.3 KB
Line 
1
2
3
4
5
6
7
8###############################################
9#OBSOLETE STUFF
10#Native checkpoint format.
11#Information needed to recreate a state is preserved
12#FIXME: Rethink and maybe use netcdf format
13def cpt_variable_writer(filename, t, v0, v1, v2):
14    """Store all conserved quantities to file
15    """
16
17    M, N = v0.shape
18
19    FN = create_filename(filename, 'cpt', M, t)
20    #print 'Writing to %s' %FN
21
22    fid = open(FN, 'w')
23    for i in range(M):
24        for j in range(N):
25            fid.write('%.16e ' %v0[i,j])
26        for j in range(N):
27            fid.write('%.16e ' %v1[i,j])
28        for j in range(N):
29            fid.write('%.16e ' %v2[i,j])
30
31        fid.write('\n')
32    fid.close()
33
34
35def cpt_variable_reader(filename, t, v0, v1, v2):
36    """Store all conserved quantities to file
37    """
38
39    M, N = v0.shape
40
41    FN = create_filename(filename, 'cpt', M, t)
42    #print 'Reading from %s' %FN
43
44    fid = open(FN)
45
46
47    for i in range(M):
48        values = fid.readline().split() #Get one line
49
50        for j in range(N):
51            v0[i,j] = float(values[j])
52            v1[i,j] = float(values[3+j])
53            v2[i,j] = float(values[6+j])
54
55    fid.close()
56
57def cpt_constant_writer(filename, X0, X1, X2, v0, v1, v2):
58    """Writes x,y,z,z,z coordinates of triangles constituting the bed
59    elevation.
60    FIXME: Not in use pt
61    """
62
63    M, N = v0.shape
64
65
66    print X0
67    import sys; sys.exit()
68    FN = create_filename(filename, 'cpt', M)
69    print 'Writing to %s' %FN
70
71    fid = open(FN, 'w')
72    for i in range(M):
73        for j in range(2):
74            fid.write('%.16e ' %X0[i,j])   #x, y
75        for j in range(N):
76            fid.write('%.16e ' %v0[i,j])       #z,z,z,
77
78        for j in range(2):
79            fid.write('%.16e ' %X1[i,j])   #x, y
80        for j in range(N):
81            fid.write('%.16e ' %v1[i,j])
82
83        for j in range(2):
84            fid.write('%.16e ' %X2[i,j])   #x, y
85        for j in range(N):
86            fid.write('%.16e ' %v2[i,j])
87
88        fid.write('\n')
89    fid.close()
90
91
92
93#Function for storing out to e.g. visualisation
94#FIXME: Do we want this?
95#FIXME: Not done yet for this version
96def dat_constant_writer(filename, X0, X1, X2, v0, v1, v2):
97    """Writes x,y,z coordinates of triangles constituting the bed elevation.
98    """
99
100    M, N = v0.shape
101
102    FN = create_filename(filename, 'dat', M)
103    #print 'Writing to %s' %FN
104
105    fid = open(FN, 'w')
106    for i in range(M):
107        for j in range(2):
108            fid.write('%f ' %X0[i,j])   #x, y
109        fid.write('%f ' %v0[i,0])       #z
110
111        for j in range(2):
112            fid.write('%f ' %X1[i,j])   #x, y
113        fid.write('%f ' %v1[i,0])       #z
114
115        for j in range(2):
116            fid.write('%f ' %X2[i,j])   #x, y
117        fid.write('%f ' %v2[i,0])       #z
118
119        fid.write('\n')
120    fid.close()
121
122
123
124def dat_variable_writer(filename, t, v0, v1, v2):
125    """Store water height to file
126    """
127
128    M, N = v0.shape
129
130    FN = create_filename(filename, 'dat', M, t)
131    #print 'Writing to %s' %FN
132
133    fid = open(FN, 'w')
134    for i in range(M):
135        fid.write('%.4f ' %v0[i,0])
136        fid.write('%.4f ' %v1[i,0])
137        fid.write('%.4f ' %v2[i,0])
138
139        fid.write('\n')
140    fid.close()
141
142
143def read_sww(filename):
144    """Read sww Net CDF file containing Shallow Water Wave simulation
145
146    The integer array volumes is of shape Nx3 where N is the number of
147    triangles in the mesh.
148
149    Each entry in volumes is an index into the x,y arrays (the location).
150
151    Quantities stage, elevation, xmomentum and ymomentum are all in arrays of dimensions
152    number_of_timesteps, number_of_points.
153
154    The momentum is not always stored.
155
156    """
157    from Scientific.IO.NetCDF import NetCDFFile
158    print 'Reading from ', filename
159    fid = NetCDFFile(filename, 'r')    #Open existing file for read
160#latitude, longitude
161    # Get the variables as Numeric arrays
162    x = fid.variables['x']             #x-coordinates of vertices
163    y = fid.variables['y']             #y-coordinates of vertices
164    z = fid.variables['elevation']     #Elevation
165    time = fid.variables['time']       #Timesteps
166    stage = fid.variables['stage']     #Water level
167    #xmomentum = fid.variables['xmomentum']   #Momentum in the x-direction
168    #ymomentum = fid.variables['ymomentum']   #Momentum in the y-direction
169
170    volumes = fid.variables['volumes'] #Connectivity
171
172    #FIXME (Ole): What is this?
173    #             Why isn't anything returned?
174    #             Where's the unit test?
175
176
177
178def sww2asc_obsolete(basename_in, basename_out = None,
179            quantity = None,
180            timestep = None,
181            reduction = None,
182            cellsize = 10,
183            verbose = False,
184            origin = None):
185    """Read SWW file and convert to Digitial Elevation model format (.asc)
186
187    Example:
188
189    ncols         3121
190    nrows         1800
191    xllcorner     722000
192    yllcorner     5893000
193    cellsize      25
194    NODATA_value  -9999
195    138.3698 137.4194 136.5062 135.5558 ..........
196
197    Also write accompanying file with same basename_in but extension .prj
198    used to fix the UTM zone, datum, false northings and eastings.
199
200    The prj format is assumed to be as
201
202    Projection    UTM
203    Zone          56
204    Datum         WGS84
205    Zunits        NO
206    Units         METERS
207    Spheroid      WGS84
208    Xshift        0.0000000000
209    Yshift        10000000.0000000000
210    Parameters
211
212
213    if quantity is given, out values from quantity otherwise default to
214    elevation
215
216    if timestep (an index) is given, output quantity at that timestep
217
218    if reduction is given use that to reduce quantity over all timesteps.
219
220    """
221    from Numeric import array, Float, concatenate, NewAxis, zeros,\
222         sometrue
223    from utilities.polygon import inside_polygon
224
225    #FIXME: Should be variable
226    datum = 'WGS84'
227    false_easting = 500000
228    false_northing = 10000000
229
230    if quantity is None:
231        quantity = 'elevation'
232
233    if reduction is None:
234        reduction = max
235
236    if basename_out is None:
237        basename_out = basename_in + '_%s' %quantity
238
239    swwfile = basename_in + '.sww'
240    ascfile = basename_out + '.asc'
241    prjfile = basename_out + '.prj'
242
243
244    if verbose: print 'Reading from %s' %swwfile
245    #Read sww file
246    from Scientific.IO.NetCDF import NetCDFFile
247    fid = NetCDFFile(swwfile)
248
249    #Get extent and reference
250    x = fid.variables['x'][:]
251    y = fid.variables['y'][:]
252    volumes = fid.variables['volumes'][:]
253
254    ymin = min(y); ymax = max(y)
255    xmin = min(x); xmax = max(x)
256
257    number_of_timesteps = fid.dimensions['number_of_timesteps']
258    number_of_points = fid.dimensions['number_of_points']
259    if origin is None:
260
261        #Get geo_reference
262        #sww files don't have to have a geo_ref
263        try:
264            geo_reference = Geo_reference(NetCDFObject=fid)
265        except AttributeError, e:
266            geo_reference = Geo_reference() #Default georef object
267
268        xllcorner = geo_reference.get_xllcorner()
269        yllcorner = geo_reference.get_yllcorner()
270        zone = geo_reference.get_zone()
271    else:
272        zone = origin[0]
273        xllcorner = origin[1]
274        yllcorner = origin[2]
275
276
277    #Get quantity and reduce if applicable
278    if verbose: print 'Reading quantity %s' %quantity
279
280    if quantity.lower() == 'depth':
281        q = fid.variables['stage'][:] - fid.variables['elevation'][:]
282    else:
283        q = fid.variables[quantity][:]
284
285
286    if len(q.shape) == 2:
287        if verbose: print 'Reducing quantity %s' %quantity
288        q_reduced = zeros( number_of_points, Float )
289
290        for k in range(number_of_points):
291            q_reduced[k] = reduction( q[:,k] )
292
293        q = q_reduced
294
295    #Now q has dimension: number_of_points
296
297    #Create grid and update xll/yll corner and x,y
298    if verbose: print 'Creating grid'
299    ncols = int((xmax-xmin)/cellsize)+1
300    nrows = int((ymax-ymin)/cellsize)+1
301
302    newxllcorner = xmin+xllcorner
303    newyllcorner = ymin+yllcorner
304
305    x = x+xllcorner-newxllcorner
306    y = y+yllcorner-newyllcorner
307
308    vertex_points = concatenate ((x[:, NewAxis] ,y[:, NewAxis]), axis = 1)
309    assert len(vertex_points.shape) == 2
310
311
312    from Numeric import zeros, Float
313    grid_points = zeros ( (ncols*nrows, 2), Float )
314
315
316    for i in xrange(nrows):
317        yg = i*cellsize
318        for j in xrange(ncols):
319            xg = j*cellsize
320            k = i*ncols + j
321
322            grid_points[k,0] = xg
323            grid_points[k,1] = yg
324
325    #Interpolate
326    from least_squares import Interpolation
327
328
329    #FIXME: This should be done with precrop = True, otherwise it'll
330    #take forever. With expand_search  set to False, some grid points might
331    #miss out....
332
333    interp = Interpolation(vertex_points, volumes, grid_points, alpha=0.0,
334                           precrop = False, expand_search = True,
335                           verbose = verbose)
336
337    #Interpolate using quantity values
338    if verbose: print 'Interpolating'
339    grid_values = interp.interpolate(q).flat
340
341    #Write
342    #Write prj file
343    if verbose: print 'Writing %s' %prjfile
344    prjid = open(prjfile, 'w')
345    prjid.write('Projection    %s\n' %'UTM')
346    prjid.write('Zone          %d\n' %zone)
347    prjid.write('Datum         %s\n' %datum)
348    prjid.write('Zunits        NO\n')
349    prjid.write('Units         METERS\n')
350    prjid.write('Spheroid      %s\n' %datum)
351    prjid.write('Xshift        %d\n' %false_easting)
352    prjid.write('Yshift        %d\n' %false_northing)
353    prjid.write('Parameters\n')
354    prjid.close()
355
356
357
358    if verbose: print 'Writing %s' %ascfile
359    NODATA_value = -9999
360
361    ascid = open(ascfile, 'w')
362
363    ascid.write('ncols         %d\n' %ncols)
364    ascid.write('nrows         %d\n' %nrows)
365    ascid.write('xllcorner     %d\n' %newxllcorner)
366    ascid.write('yllcorner     %d\n' %newyllcorner)
367    ascid.write('cellsize      %f\n' %cellsize)
368    ascid.write('NODATA_value  %d\n' %NODATA_value)
369
370
371    #Get bounding polygon from mesh
372    P = interp.mesh.get_boundary_polygon()
373    inside_indices = inside_polygon(grid_points, P)
374
375    for i in range(nrows):
376        if verbose and i%((nrows+10)/10)==0:
377            print 'Doing row %d of %d' %(i, nrows)
378
379        for j in range(ncols):
380            index = (nrows-i-1)*ncols+j
381
382            if sometrue(inside_indices == index):
383                ascid.write('%f ' %grid_values[index])
384            else:
385                ascid.write('%d ' %NODATA_value)
386
387        ascid.write('\n')
388
389    #Close
390    ascid.close()
391    fid.close()
392
393#********************
394#*** END OF OBSOLETE FUNCTIONS
395#***************
Note: See TracBrowser for help on using the repository browser.