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

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

Fixed up failing tests, updated user guide with new API (first few chapters only).

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, urs_ungridded2sww
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.