source: branches/numpy/obsolete_code/data_manager_obsolete_stuff.py @ 8353

Last change on this file since 8353 was 3761, checked in by ole, 18 years ago

Removed old tests in data_manager as per ticket:86 - they were old versions of dem2pts tests that had been superseded. All good.
File size: 21.5 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 anuga.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
395
396    def NOT_test_dem2pts(self):
397        """Test conversion from dem in ascii format to native NetCDF xya format
398        """
399
400        import time, os
401        from Numeric import array, zeros, allclose, Float, concatenate
402        from Scientific.IO.NetCDF import NetCDFFile
403
404        #Write test asc file
405        root = 'demtest'
406
407        filename = root+'.asc'
408        fid = open(filename, 'w')
409        fid.write("""ncols         5
410nrows         6
411xllcorner     2000.5
412yllcorner     3000.5
413cellsize      25
414NODATA_value  -9999
415""")
416        #Create linear function
417
418        ref_points = []
419        ref_elevation = []
420        for i in range(6):
421            y = (6-i)*25.0
422            for j in range(5):
423                x = j*25.0
424                z = x+2*y
425
426                ref_points.append( [x,y] )
427                ref_elevation.append(z)
428                fid.write('%f ' %z)
429            fid.write('\n')
430
431        fid.close()
432
433        #Write prj file with metadata
434        metafilename = root+'.prj'
435        fid = open(metafilename, 'w')
436
437
438        fid.write("""Projection UTM
439Zone 56
440Datum WGS84
441Zunits NO
442Units METERS
443Spheroid WGS84
444Xshift 0.0000000000
445Yshift 10000000.0000000000
446Parameters
447""")
448        fid.close()
449
450        #Convert to NetCDF pts
451        convert_dem_from_ascii2netcdf(root)
452        dem2pts(root)
453
454        #Check contents
455        #Get NetCDF
456        fid = NetCDFFile(root+'.pts', 'r')
457
458        # Get the variables
459        #print fid.variables.keys()
460        points = fid.variables['points']
461        elevation = fid.variables['elevation']
462
463        #Check values
464
465        #print points[:]
466        #print ref_points
467        assert allclose(points, ref_points)
468
469        #print attributes[:]
470        #print ref_elevation
471        assert allclose(elevation, ref_elevation)
472
473        #Cleanup
474        fid.close()
475
476
477        os.remove(root + '.pts')
478        os.remove(root + '.dem')
479        os.remove(root + '.asc')
480        os.remove(root + '.prj')
481
482
483
484    def NOT_test_dem2pts_bounding_box(self):
485        """Test conversion from dem in ascii format to native NetCDF xya format
486        """
487
488        import time, os
489        from Numeric import array, zeros, allclose, Float, concatenate
490        from Scientific.IO.NetCDF import NetCDFFile
491
492        #Write test asc file
493        root = 'demtest'
494
495        filename = root+'.asc'
496        fid = open(filename, 'w')
497        fid.write("""ncols         5
498nrows         6
499xllcorner     2000.5
500yllcorner     3000.5
501cellsize      25
502NODATA_value  -9999
503""")
504        #Create linear function
505
506        ref_points = []
507        ref_elevation = []
508        for i in range(6):
509            y = (6-i)*25.0
510            for j in range(5):
511                x = j*25.0
512                z = x+2*y
513
514                ref_points.append( [x,y] )
515                ref_elevation.append(z)
516                fid.write('%f ' %z)
517            fid.write('\n')
518
519        fid.close()
520
521        #Write prj file with metadata
522        metafilename = root+'.prj'
523        fid = open(metafilename, 'w')
524
525
526        fid.write("""Projection UTM
527Zone 56
528Datum WGS84
529Zunits NO
530Units METERS
531Spheroid WGS84
532Xshift 0.0000000000
533Yshift 10000000.0000000000
534Parameters
535""")
536        fid.close()
537
538        #Convert to NetCDF pts
539        convert_dem_from_ascii2netcdf(root)
540        dem2pts(root, easting_min=2010.0, easting_max=2110.0,
541                northing_min=3035.0, northing_max=3125.5)
542
543        #Check contents
544        #Get NetCDF
545        fid = NetCDFFile(root+'.pts', 'r')
546
547        # Get the variables
548        #print fid.variables.keys()
549        points = fid.variables['points']
550        elevation = fid.variables['elevation']
551
552        #Check values
553        assert fid.xllcorner[0] == 2010.0
554        assert fid.yllcorner[0] == 3035.0
555
556        #create new reference points
557        ref_points = []
558        ref_elevation = []
559        for i in range(4):
560            y = (4-i)*25.0 + 25.0
561            y_new = y + 3000.5 - 3035.0
562            for j in range(4):
563                x = j*25.0 + 25.0
564                x_new = x + 2000.5 - 2010.0
565                z = x+2*y
566
567                ref_points.append( [x_new,y_new] )
568                ref_elevation.append(z)
569
570        #print points[:]
571        #print ref_points
572        assert allclose(points, ref_points)
573
574        #print attributes[:]
575        #print ref_elevation
576        assert allclose(elevation, ref_elevation)
577
578        #Cleanup
579        fid.close()
580
581
582        os.remove(root + '.pts')
583        os.remove(root + '.dem')
584        os.remove(root + '.asc')
585        os.remove(root + '.prj')
586
587
588
589    def NOT_test_dem2pts_remove_Nullvalues(self):
590        """Test conversion from dem in ascii format to native NetCDF xya format
591        """
592
593        import time, os
594        from Numeric import array, zeros, allclose, Float, concatenate
595        from Scientific.IO.NetCDF import NetCDFFile
596
597        #Write test asc file
598        root = 'demtest'
599
600        filename = root+'.asc'
601        fid = open(filename, 'w')
602        fid.write("""ncols         5
603nrows         6
604xllcorner     2000.5
605yllcorner     3000.5
606cellsize      25
607NODATA_value  -9999
608""")
609        #Create linear function
610        # ref_ will write all the values
611        # new_ref_ will write the values except for NODATA_values
612        ref_points = []
613        ref_elevation = []
614        new_ref_pts = []
615        new_ref_elev = []
616        NODATA_value = -9999
617        for i in range(6):
618            y = (6-i)*25.0
619            for j in range(5):
620                x = j*25.0
621                z = x+2*y
622                if j == 4: z = NODATA_value # column
623                if i == 2 and j == 2: z = NODATA_value # random
624                if i == 5 and j == 1: z = NODATA_value
625                if i == 1: z = NODATA_value # row
626                if i == 3 and j == 1: z = NODATA_value # two pts/row
627                if i == 3 and j == 3: z = NODATA_value
628
629
630                if z <> NODATA_value:
631                    new_ref_elev.append(z)
632                    new_ref_pts.append( [x,y] )
633
634                ref_points.append( [x,y] )
635                ref_elevation.append(z)
636
637                fid.write('%f ' %z)
638            fid.write('\n')
639
640        fid.close()
641
642
643        #Write prj file with metadata
644        metafilename = root+'.prj'
645        fid = open(metafilename, 'w')
646
647
648        fid.write("""Projection UTM
649Zone 56
650Datum WGS84
651Zunits NO
652Units METERS
653Spheroid WGS84
654Xshift 0.0000000000
655Yshift 10000000.0000000000
656Parameters
657""")
658        fid.close()
659
660        #Convert to NetCDF pts
661        convert_dem_from_ascii2netcdf(root)
662        dem2pts(root)
663
664        #Check contents
665        #Get NetCDF
666        fid = NetCDFFile(root+'.pts', 'r')
667
668        # Get the variables
669        #print fid.variables.keys()
670        points = fid.variables['points']
671        elevation = fid.variables['elevation']
672
673        #Check values
674        #print 'points', points[:]
675        assert len(points) == len(new_ref_pts), 'length of returned points not correct'
676        assert allclose(points, new_ref_pts), 'points do not align'
677
678        #print 'elevation', elevation[:]
679        assert len(elevation) == len(new_ref_elev), 'length of returned elevation not correct'
680        assert allclose(elevation, new_ref_elev), 'elevations do not align'
681
682        #Cleanup
683        fid.close()
684
685
686        os.remove(root + '.pts')
687        os.remove(root + '.dem')
688        os.remove(root + '.asc')
689        os.remove(root + '.prj')
690
691    def NOT_test_dem2pts_bounding_box_Nullvalues(self):
692        """Test conversion from dem in ascii format to native NetCDF xya format
693        """
694
695        import time, os
696        from Numeric import array, zeros, allclose, Float, concatenate
697        from Scientific.IO.NetCDF import NetCDFFile
698
699        #Write test asc file
700        root = 'demtest'
701
702        filename = root+'.asc'
703        fid = open(filename, 'w')
704        fid.write("""ncols         5
705nrows         6
706xllcorner     2000.5
707yllcorner     3000.5
708cellsize      25
709NODATA_value  -9999
710""")
711        #Create linear function
712
713        ref_points = []
714        ref_elevation = []
715        new_ref_pts1 = []
716        new_ref_elev1 = []
717        NODATA_value = -9999
718        for i in range(6):
719            y = (6-i)*25.0
720            for j in range(5):
721                x = j*25.0
722                z = x+2*y
723                if j == 4: z = NODATA_value # column
724                if i == 2 and j == 2: z = NODATA_value # random
725                if i == 5 and j == 1: z = NODATA_value
726                if i == 1: z = NODATA_value # row
727                if i == 3 and j == 1: z = NODATA_value # two pts/row
728                if i == 3 and j == 3: z = NODATA_value
729
730                if z <> NODATA_value:
731                    new_ref_elev1.append(z)
732                    new_ref_pts1.append( [x,y] )
733
734                ref_points.append( [x,y] )
735                ref_elevation.append(z)
736                fid.write('%f ' %z)
737            fid.write('\n')
738
739        fid.close()
740
741        #Write prj file with metadata
742        metafilename = root+'.prj'
743        fid = open(metafilename, 'w')
744
745
746        fid.write("""Projection UTM
747Zone 56
748Datum WGS84
749Zunits NO
750Units METERS
751Spheroid WGS84
752Xshift 0.0000000000
753Yshift 10000000.0000000000
754Parameters
755""")
756        fid.close()
757
758        #Convert to NetCDF pts
759        convert_dem_from_ascii2netcdf(root)
760        dem2pts(root, easting_min=2010.0, easting_max=2110.0,
761                northing_min=3035.0, northing_max=3125.5)
762
763        #Check contents
764        #Get NetCDF
765        fid = NetCDFFile(root+'.pts', 'r')
766
767        # Get the variables
768        #print fid.variables.keys()
769        points = fid.variables['points']
770        elevation = fid.variables['elevation']
771
772        #Check values
773        assert fid.xllcorner[0] == 2010.0
774        assert fid.yllcorner[0] == 3035.0
775
776        #create new reference points
777        ref_points = []
778        ref_elevation = []
779        new_ref_pts2 = []
780        new_ref_elev2 = []
781        for i in range(4):
782            y = (4-i)*25.0 + 25.0
783            y_new = y + 3000.5 - 3035.0
784            for j in range(4):
785                x = j*25.0 + 25.0
786                x_new = x + 2000.5 - 2010.0
787                z = x+2*y
788
789                if j == 3: z = NODATA_value # column
790                if i == 1 and j == 1: z = NODATA_value # random
791                if i == 4 and j == 0: z = NODATA_value
792                if i == 0: z = NODATA_value # row
793                if i == 2 and j == 0: z = NODATA_value # two pts/row
794                if i == 2 and j == 2: z = NODATA_value
795
796                if z <> NODATA_value:
797                    new_ref_elev2.append(z)
798                    new_ref_pts2.append( [x_new,y_new] )
799
800
801                ref_points.append( [x_new,y_new] )
802                ref_elevation.append(z)
803
804        #print points[:]
805        #print ref_points
806        #assert allclose(points, ref_points)
807
808        #print attributes[:]
809        #print ref_elevation
810        #assert allclose(elevation, ref_elevation)
811
812
813        assert len(points) == len(new_ref_pts2), 'length of returned points not correct'
814        assert allclose(points, new_ref_pts2), 'points do not align'
815
816        #print 'elevation', elevation[:]
817        assert len(elevation) == len(new_ref_elev2), 'length of returned elevation not correct'
818        assert allclose(elevation, new_ref_elev2), 'elevations do not align'
819        #Cleanup
820        fid.close()
821
822
823        os.remove(root + '.pts')
824        os.remove(root + '.dem')
825        os.remove(root + '.asc')
826        os.remove(root + '.prj')
827
828#********************
829#*** END OF OBSOLETE FUNCTIONS
830#***************
Note: See TracBrowser for help on using the repository browser.