source: trunk/anuga_core/source/anuga/file_conversion/test_urs2sww.py @ 7800

Last change on this file since 7800 was 7800, checked in by hudson, 12 years ago

urs2sww has an extra urs_ungridded2sww function.

File size: 23.8 KB
Line 
1
2# External modules
3from Scientific.IO.NetCDF import NetCDFFile
4import sys
5import unittest
6import numpy as num
7import copy
8import os
9
10# ANUGA modules
11from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a, \
12                            netcdf_float
13
14from anuga.coordinate_transforms.geo_reference import Geo_reference, \
15     write_NetCDF_georeference
16from anuga.coordinate_transforms.redfearn import redfearn     
17
18from urs2sww import urs2sww
19import urs
20
21from anuga.file.mux import WAVEHEIGHT_MUX_LABEL, EAST_VELOCITY_LABEL, \
22                            NORTH_VELOCITY_LABEL
23
24from anuga.geospatial_data.geospatial_data import ensure_absolute
25from anuga.utilities.numerical_tools import ensure_numeric 
26
27# use helper methods from other unit test
28from anuga.file.test_mux import Test_Mux
29
30
31
32class Test_Dem2Pts(Test_Mux):
33    """ A suite of tests to test urs2sww file conversion functions.
34        These tests are quite coarse-grained: converting a file
35        and checking that its headers and some of its contents
36        are correct.
37    """ 
38
39    def setUp(self):
40        self.verbose = False # change this to output more debug info
41
42 
43    def test_urs_ungridded2swwIII (self):
44       
45        #Zone:   50   
46        #Easting:  240992.578  Northing: 7620442.472
47        #Latitude:   -21  30 ' 0.00000 ''  Longitude: 114  30 ' 0.00000 ''
48        lat_long = [[-21.5,114.5],[-21,114.5],[-21,115]]
49        time_step_count = 2
50        time_step = 400
51        tide = 9000000
52        base_name, files = self.write_mux(lat_long,
53                                          time_step_count, time_step)
54        urs_ungridded2sww(base_name, mean_stage=tide, origin =(50,23432,4343))
55       
56        # now I want to check the sww file ...
57        sww_file = base_name + '.sww'
58       
59        #Let's interigate the sww file
60        # Note, the sww info is not gridded.  It is point data.
61        fid = NetCDFFile(sww_file)
62       
63        # Make x and y absolute
64        x = fid.variables['x'][:]
65        y = fid.variables['y'][:]
66        geo_reference = Geo_reference(NetCDFObject=fid)
67        points = geo_reference.get_absolute(map(None, x, y))
68        points = ensure_numeric(points)
69        x = points[:,0]
70        y = points[:,1]
71       
72        #Check that first coordinate is correctly represented       
73        #Work out the UTM coordinates for first point
74        zone, e, n = redfearn(lat_long[0][0], lat_long[0][1]) 
75        assert num.allclose([x[0],y[0]], [e,n])
76
77        #Check the time vector
78        times = fid.variables['time'][:]
79       
80        times_actual = []
81        for i in range(time_step_count):
82            times_actual.append(time_step * i)
83       
84        assert num.allclose(ensure_numeric(times),
85                            ensure_numeric(times_actual))
86       
87        #Check first value
88        stage = fid.variables['stage'][:]
89        xmomentum = fid.variables['xmomentum'][:]
90        ymomentum = fid.variables['ymomentum'][:]
91        elevation = fid.variables['elevation'][:]
92        assert num.allclose(stage[0,0], e +tide)  #Meters
93
94        #Check the momentums - ua
95        #momentum = velocity*(stage-elevation)
96        # elevation = - depth
97        #momentum = velocity_ua *(stage+depth)
98        # = n*(e+tide+n) based on how I'm writing these files
99        #
100        answer_x = n*(e+tide+n)
101        actual_x = xmomentum[0,0]
102        #print "answer_x",answer_x
103        #print "actual_x",actual_x
104        assert num.allclose(answer_x, actual_x)  #Meters
105       
106        #Check the momentums - va
107        #momentum = velocity*(stage-elevation)
108        # elevation = - depth
109        #momentum = velocity_va *(stage+depth)
110        # = e*(e+tide+n) based on how I'm writing these files
111        #
112        answer_y = -1*e*(e+tide+n)
113        actual_y = ymomentum[0,0]
114        #print "answer_y",answer_y
115        #print "actual_y",actual_y
116        assert num.allclose(answer_y, actual_y)  #Meters
117
118        # check the stage values, first time step.
119        # These arrays are equal since the Easting values were used as
120        # the stage
121        assert num.allclose(stage[0], x +tide)  #Meters
122        # check the elevation values.
123        # -ve since urs measures depth, sww meshers height,
124        # these arrays are equal since the northing values were used as
125        # the elevation
126        assert num.allclose(-elevation, y)  #Meters
127       
128        fid.close()
129        self.delete_mux(files)
130        os.remove(sww_file)
131
132       
133    def test_urs_ungridded_hole (self):
134       
135        #Zone:   50   
136        #Easting:  240992.578  Northing: 7620442.472
137        #Latitude:   -21  30 ' 0.00000 ''  Longitude: 114  30 ' 0.00000 ''
138        lat_long = [[-20.5, 114.5],
139                    [-20.6, 114.6],
140                    [-20.5, 115.],
141                    [-20.6, 115.],
142                    [-20.5, 115.5],
143                    [-20.6, 115.4],
144                   
145                    [-21., 114.5],
146                    [-21., 114.6],
147                    [-21., 115.5],
148                    [-21., 115.4],
149                   
150                    [-21.5, 114.5],
151                    [-21.4, 114.6],
152                    [-21.5, 115.],
153                    [-21.4, 115.],
154                    [-21.5, 115.5],
155                    [-21.4, 115.4]
156                    ]
157        time_step_count = 6
158        time_step = 100
159        tide = 9000000
160        base_name, files = self.write_mux(lat_long,
161                                          time_step_count, time_step)
162        #Easting:  292110.784  Northing: 7676551.710
163        #Latitude:   -21  0 ' 0.00000 ''  Longitude: 115  0 ' 0.00000 ''
164
165        'hole_points_UTM(base_name, mean_stage=-240992.0,
166                          hole_points_UTM=[ 292110.784, 7676551.710 ])
167       
168        # now I want to check the sww file ...
169        sww_file = base_name + '.sww'
170       
171        #Let's interigate the sww file
172        # Note, the sww info is not gridded.  It is point data.
173        fid = NetCDFFile(sww_file)
174       
175        number_of_volumes = fid.variables['volumes']
176        #print "number_of_volumes",len(number_of_volumes)
177        assert num.allclose(16, len(number_of_volumes))
178       
179        fid.close()
180        self.delete_mux(files)
181        #print "sww_file", sww_file
182        os.remove(sww_file)
183       
184    def test_urs_ungridded_holeII(self):
185
186        # Check that if using a hole that returns no triangles,
187        # urs_ungridded2sww removes the hole label.
188       
189        lat_long = [[-20.5, 114.5],
190                    [-20.6, 114.6],
191                    [-20.5, 115.5],
192                    [-20.6, 115.4],
193                   
194                   
195                    [-21.5, 114.5],
196                    [-21.4, 114.6],
197                    [-21.5, 115.5],
198                    [-21.4, 115.4]
199                    ]
200        time_step_count = 6
201        time_step = 100
202        tide = 9000000
203        base_name, files = self.write_mux(lat_long,
204                                          time_step_count, time_step)
205        #Easting:  292110.784  Northing: 7676551.710
206        #Latitude:   -21  0 ' 0.00000 ''  Longitude: 115  0 ' 0.00000 ''
207
208        urs_ungridded2sww(base_name, mean_stage=-240992.0,
209                          hole_points_UTM=[ 292110.784, 7676551.710 ])
210       
211        # now I want to check the sww file ...
212        sww_file = base_name + '.sww'
213        fid = NetCDFFile(sww_file)
214       
215        volumes = fid.variables['volumes']
216        #print "number_of_volumes",len(volumes)
217
218        fid.close()
219        os.remove(sww_file)
220       
221        urs2sww(base_name, mean_stage=-240992.0)
222       
223        # now I want to check the sww file ...
224        sww_file = base_name + '.sww'
225        fid = NetCDFFile(sww_file)
226       
227        volumes_again = fid.variables['volumes']
228        #print "number_of_volumes",len(volumes_again)
229        assert num.allclose(len(volumes_again),
230                            len(volumes))
231        fid.close()
232        os.remove(sww_file)
233        self.delete_mux(files) 
234       
235    def test_urs_ungridded2sww_mint_maxt (self):
236       
237        #Zone:   50   
238        #Easting:  240992.578  Northing: 7620442.472
239        #Latitude:   -21  30 ' 0.00000 ''  Longitude: 114  30 ' 0.00000 ''
240        lat_long = [[-21.5,114.5],[-21,114.5],[-21,115]]
241        time_step_count = 6
242        time_step = 100
243        tide = 9000000
244        base_name, files = self.write_mux(lat_long,
245                                          time_step_count, time_step)
246        urs_ungridded2sww(base_name, mean_stage=tide, origin =(50,23432,4343),
247                          mint=101, maxt=500)
248       
249        # now I want to check the sww file ...
250        sww_file = base_name + '.sww'
251       
252        #Let's interigate the sww file
253        # Note, the sww info is not gridded.  It is point data.
254        fid = NetCDFFile(sww_file)
255       
256        # Make x and y absolute
257        x = fid.variables['x'][:]
258        y = fid.variables['y'][:]
259        geo_reference = Geo_reference(NetCDFObject=fid)
260        points = geo_reference.get_absolute(map(None, x, y))
261        points = ensure_numeric(points)
262        x = points[:,0]
263        y = points[:,1]
264       
265        #Check that first coordinate is correctly represented       
266        #Work out the UTM coordinates for first point
267        zone, e, n = redfearn(lat_long[0][0], lat_long[0][1]) 
268        assert num.allclose([x[0],y[0]], [e,n])
269
270        #Check the time vector
271        times = fid.variables['time'][:]
272       
273        times_actual = [0,100,200,300]
274       
275        assert num.allclose(ensure_numeric(times),
276                            ensure_numeric(times_actual))
277       
278               #Check first value
279        stage = fid.variables['stage'][:]
280        xmomentum = fid.variables['xmomentum'][:]
281        ymomentum = fid.variables['ymomentum'][:]
282        elevation = fid.variables['elevation'][:]
283        assert num.allclose(stage[0,0], e +tide)  #Meters
284
285        #Check the momentums - ua
286        #momentum = velocity*(stage-elevation)
287        # elevation = - depth
288        #momentum = velocity_ua *(stage+depth)
289        # = n*(e+tide+n) based on how I'm writing these files
290        #
291        answer_x = n*(e+tide+n)
292        actual_x = xmomentum[0,0]
293        #print "answer_x",answer_x
294        #print "actual_x",actual_x
295        assert num.allclose(answer_x, actual_x)  #Meters
296       
297        #Check the momentums - va
298        #momentum = velocity*(stage-elevation)
299        # elevation = - depth
300        #momentum = velocity_va *(stage+depth)
301        # = e*(e+tide+n) based on how I'm writing these files
302        #
303        answer_y = -1*e*(e+tide+n)
304        actual_y = ymomentum[0,0]
305        #print "answer_y",answer_y
306        #print "actual_y",actual_y
307        assert num.allclose(answer_y, actual_y)  #Meters
308
309        # check the stage values, first time step.
310        # These arrays are equal since the Easting values were used as
311        # the stage
312        assert num.allclose(stage[0], x +tide)  #Meters
313        # check the elevation values.
314        # -ve since urs measures depth, sww meshers height,
315        # these arrays are equal since the northing values were used as
316        # the elevation
317        assert num.allclose(-elevation, y)  #Meters
318       
319        fid.close()
320        self.delete_mux(files)
321        os.remove(sww_file)
322       
323    def test_urs_ungridded2sww_mint_maxtII (self):
324       
325        #Zone:   50   
326        #Easting:  240992.578  Northing: 7620442.472
327        #Latitude:   -21  30 ' 0.00000 ''  Longitude: 114  30 ' 0.00000 ''
328        lat_long = [[-21.5,114.5],[-21,114.5],[-21,115]]
329        time_step_count = 6
330        time_step = 100
331        tide = 9000000
332        base_name, files = self.write_mux(lat_long,
333                                          time_step_count, time_step)
334        urs_ungridded2sww(base_name, mean_stage=tide, origin =(50,23432,4343),
335                          mint=0, maxt=100000)
336       
337        # now I want to check the sww file ...
338        sww_file = base_name + '.sww'
339       
340        #Let's interigate the sww file
341        # Note, the sww info is not gridded.  It is point data.
342        fid = NetCDFFile(sww_file)
343       
344        # Make x and y absolute
345        geo_reference = Geo_reference(NetCDFObject=fid)
346        points = geo_reference.get_absolute(map(None, fid.variables['x'][:],
347                                                fid.variables['y'][:]))
348        points = ensure_numeric(points)
349        x = points[:,0]
350       
351        #Check the time vector
352        times = fid.variables['time'][:]
353       
354        times_actual = [0,100,200,300,400,500]
355        assert num.allclose(ensure_numeric(times),
356                            ensure_numeric(times_actual))
357       
358        #Check first value
359        stage = fid.variables['stage'][:]
360        assert num.allclose(stage[0], x +tide)
361       
362        fid.close()
363        self.delete_mux(files)
364        os.remove(sww_file)
365       
366    def test_urs_ungridded2sww_mint_maxtIII (self):
367       
368        #Zone:   50   
369        #Easting:  240992.578  Northing: 7620442.472
370        #Latitude:   -21  30 ' 0.00000 ''  Longitude: 114  30 ' 0.00000 ''
371        lat_long = [[-21.5,114.5],[-21,114.5],[-21,115]]
372        time_step_count = 6
373        time_step = 100
374        tide = 9000000
375        base_name, files = self.write_mux(lat_long,
376                                          time_step_count, time_step)
377        try:
378            urs2sww(base_name, mean_stage=tide,
379                          origin =(50,23432,4343),
380                          mint=301, maxt=399,
381                              verbose=self.verbose)
382        except: 
383            pass
384        else:
385            self.failUnless(0 ==1, 'Bad input did not throw exception error!')
386
387        self.delete_mux(files)
388       
389    def test_urs_ungridded2sww_mint_maxt_bad (self):       
390        #Zone:   50   
391        #Easting:  240992.578  Northing: 7620442.472
392        #Latitude:   -21  30 ' 0.00000 ''  Longitude: 114  30 ' 0.00000 ''
393        lat_long = [[-21.5,114.5],[-21,114.5],[-21,115]]
394        time_step_count = 6
395        time_step = 100
396        tide = 9000000
397        base_name, files = self.write_mux(lat_long,
398                                          time_step_count, time_step)
399        try:
400            urs2sww(base_name, mean_stage=tide,
401                          origin =(50,23432,4343),
402                          mint=301, maxt=301,
403                              verbose=self.verbose)
404        except: 
405            pass
406        else:
407            self.failUnless(0 ==1, 'Bad input did not throw exception error!')
408
409        self.delete_mux(files)
410
411       
412    def test_URS_points_needed_and_urs_ungridded2sww(self):
413        # This doesn't actually check anything
414       
415        ll_lat = -21.5
416        ll_long = 114.5
417        grid_spacing = 1./60.
418        lat_amount = 30
419        long_amount = 30
420        time_step_count = 2
421        time_step = 400
422        tide = -200000
423        zone = 50
424
425        boundary_polygon = [[250000,7660000],[280000,7660000],
426                             [280000,7630000],[250000,7630000]]
427        geo=urs.calculate_boundary_points(boundary_polygon, zone,
428                              ll_lat, ll_long, grid_spacing, 
429                              lat_amount, long_amount,
430                              verbose=self.verbose)
431        lat_long = geo.get_data_points(as_lat_long=True)
432        base_name, files = self.write_mux(lat_long,
433                                          time_step_count, time_step)
434        urs_ungridded2sww(base_name, mean_stage=tide,
435                          verbose=self.verbose)
436        self.delete_mux(files)
437        os.remove( base_name + '.sww')
438   
439    def cache_test_URS_points_needed_and_urs_ungridded2sww(self):
440       
441        ll_lat = -21.5
442        ll_long = 114.5
443        grid_spacing = 1./60.
444        lat_amount = 30
445        long_amount = 30
446        time_step_count = 2
447        time_step = 400
448        tide = -200000
449        zone = 50
450
451        boundary_polygon = [[250000,7660000],[270000,7650000],
452                             [280000,7630000],[250000,7630000]]
453        geo=URS_points_needed(boundary_polygon, zone,
454                              ll_lat, ll_long, grid_spacing, 
455                              lat_amount, long_amount, use_cache=True,
456                              verbose=True)
457       
458    def visual_test_URS_points_needed_and_urs_ungridded2sww(self):
459       
460        ll_lat = -21.5
461        ll_long = 114.5
462        grid_spacing = 1./60.
463        lat_amount = 30
464        long_amount = 30
465        time_step_count = 2
466        time_step = 400
467        tide = -200000
468        zone = 50
469
470        boundary_polygon = [[250000,7660000],[270000,7650000],
471                             [280000,7630000],[250000,7630000]]
472        geo=URS_points_needed(boundary_polygon, zone,
473                              ll_lat, ll_long, grid_spacing, 
474                              lat_amount, long_amount)
475        lat_long = geo.get_data_points(as_lat_long=True)
476        base_name, files = self.write_mux(lat_long,
477                                          time_step_count, time_step)
478        urs2sww(base_name, mean_stage=tide)
479        self.delete_mux(files)
480        os.remove( base_name + '.sww')
481        # extend this so it interpolates onto the boundary.
482        # have it fail if there is NaN
483
484
485 
486    def test_urs_ungridded2swwII (self):
487       
488        #Zone:   50   
489        #Easting:  240992.578  Northing: 7620442.472
490        #Latitude:   -21  30 ' 0.00000 ''  Longitude: 114  30 ' 0.00000 ''
491        lat_long = [[-21.5,114.5],[-21,114.5],[-21,115]]
492        time_step_count = 2
493        time_step = 400
494        tide = 9000000
495        geo_reference = Geo_reference(50, 3434543,34534543)
496        base_name, files = self.write_mux(lat_long,
497                                          time_step_count, time_step)
498        urs_ungridded2sww(base_name, mean_stage=tide, origin = geo_reference)
499       
500        # now I want to check the sww file ...
501        sww_file = base_name + '.sww'
502       
503        #Let's interigate the sww file
504        # Note, the sww info is not gridded.  It is point data.
505        fid = NetCDFFile(sww_file)
506       
507        # Make x and y absolute
508        x = fid.variables['x'][:]
509        y = fid.variables['y'][:]
510        geo_reference = Geo_reference(NetCDFObject=fid)
511        points = geo_reference.get_absolute(map(None, x, y))
512        points = ensure_numeric(points)
513        x = points[:,0]
514        y = points[:,1]
515       
516        #Check that first coordinate is correctly represented       
517        #Work out the UTM coordinates for first point
518        zone, e, n = redfearn(lat_long[0][0], lat_long[0][1]) 
519        assert num.allclose([x[0],y[0]], [e,n])
520
521        #Check the time vector
522        times = fid.variables['time'][:]
523       
524        times_actual = []
525        for i in range(time_step_count):
526            times_actual.append(time_step * i)
527       
528        assert num.allclose(ensure_numeric(times),
529                            ensure_numeric(times_actual))
530       
531        #Check first value
532        stage = fid.variables['stage'][:]
533        xmomentum = fid.variables['xmomentum'][:]
534        ymomentum = fid.variables['ymomentum'][:]
535        elevation = fid.variables['elevation'][:]
536        assert num.allclose(stage[0,0], e +tide)  #Meters
537
538        #Check the momentums - ua
539        #momentum = velocity*(stage-elevation)
540        # elevation = - depth
541        #momentum = velocity_ua *(stage+depth)
542        # = n*(e+tide+n) based on how I'm writing these files
543        #
544        answer_x = n*(e+tide+n)
545        actual_x = xmomentum[0,0]
546        #print "answer_x",answer_x
547        #print "actual_x",actual_x
548        assert num.allclose(answer_x, actual_x)  #Meters
549       
550        #Check the momentums - va
551        #momentum = velocity*(stage-elevation)
552        # elevation = - depth
553        #momentum = velocity_va *(stage+depth)
554        # = e*(e+tide+n) based on how I'm writing these files
555        #
556        answer_y = -1*e*(e+tide+n)
557        actual_y = ymomentum[0,0]
558        #print "answer_y",answer_y
559        #print "actual_y",actual_y
560        assert num.allclose(answer_y, actual_y)  #Meters
561
562        # check the stage values, first time step.
563        # These arrays are equal since the Easting values were used as
564        # the stage
565        assert num.allclose(stage[0], x +tide)  #Meters
566        # check the elevation values.
567        # -ve since urs measures depth, sww meshers height,
568        # these arrays are equal since the northing values were used as
569        # the elevation
570        assert num.allclose(-elevation, y)  #Meters
571       
572        fid.close()
573        self.delete_mux(files)
574        os.remove(sww_file)
575
576
577    def test_urs_ungridded2sww (self):
578       
579        #Zone:   50   
580        #Easting:  240992.578  Northing: 7620442.472
581        #Latitude:   -21  30 ' 0.00000 ''  Longitude: 114  30 ' 0.00000 ''
582        lat_long = [[-21.5,114.5],[-21,114.5],[-21,115]]
583        time_step_count = 2
584        time_step = 400
585        tide = 9000000
586        base_name, files = self.write_mux(lat_long,
587                                          time_step_count, time_step)
588        urs_ungridded2sww(base_name, mean_stage=tide,
589                          verbose=self.verbose)
590       
591        # now I want to check the sww file ...
592        sww_file = base_name + '.sww'
593       
594        #Let's interigate the sww file
595        # Note, the sww info is not gridded.  It is point data.
596        fid = NetCDFFile(sww_file)
597       
598        # Make x and y absolute
599        x = fid.variables['x'][:]
600        y = fid.variables['y'][:]
601        geo_reference = Geo_reference(NetCDFObject=fid)
602        points = geo_reference.get_absolute(map(None, x, y))
603        points = ensure_numeric(points)
604        x = points[:,0]
605        y = points[:,1]
606       
607        #Check that first coordinate is correctly represented       
608        #Work out the UTM coordinates for first point
609        zone, e, n = redfearn(lat_long[0][0], lat_long[0][1]) 
610        assert num.allclose([x[0],y[0]], [e,n])
611
612        #Check the time vector
613        times = fid.variables['time'][:]
614       
615        times_actual = []
616        for i in range(time_step_count):
617            times_actual.append(time_step * i)
618       
619        assert num.allclose(ensure_numeric(times),
620                            ensure_numeric(times_actual))
621       
622        #Check first value
623        stage = fid.variables['stage'][:]
624        xmomentum = fid.variables['xmomentum'][:]
625        ymomentum = fid.variables['ymomentum'][:]
626        elevation = fid.variables['elevation'][:]
627        assert num.allclose(stage[0,0], e +tide)  #Meters
628
629
630        #Check the momentums - ua
631        #momentum = velocity*(stage-elevation)
632        # elevation = - depth
633        #momentum = velocity_ua *(stage+depth)
634        # = n*(e+tide+n) based on how I'm writing these files
635        #
636        answer_x = n*(e+tide+n)
637        actual_x = xmomentum[0,0]
638        #print "answer_x",answer_x
639        #print "actual_x",actual_x
640        assert num.allclose(answer_x, actual_x)  #Meters
641       
642        #Check the momentums - va
643        #momentum = velocity*(stage-elevation)
644        # elevation = - depth
645        #momentum = velocity_va *(stage+depth)
646        # = e*(e+tide+n) based on how I'm writing these files
647        #
648        answer_y = -1*e*(e+tide+n)
649        actual_y = ymomentum[0,0]
650        #print "answer_y",answer_y
651        #print "actual_y",actual_y
652        assert num.allclose(answer_y, actual_y)  #Meters
653
654        # check the stage values, first time step.
655        # These arrays are equal since the Easting values were used as
656        # the stage
657        assert num.allclose(stage[0], x +tide)  #Meters
658        # check the elevation values.
659        # -ve since urs measures depth, sww meshers height,
660        # these arrays are equal since the northing values were used as
661        # the elevation
662        assert num.allclose(-elevation, y)  #Meters
663       
664        fid.close()
665        self.delete_mux(files)
666        os.remove(sww_file)
667
668
669#-------------------------------------------------------------
670
671if __name__ == "__main__":
672    suite = unittest.makeSuite(Test_Dem2Pts,'test')
673    runner = unittest.TextTestRunner()
674    runner.run(suite)
Note: See TracBrowser for help on using the repository browser.