source: trunk/anuga_core/source/anuga/file/test_urs.py @ 8813

Last change on this file since 8813 was 8249, checked in by steve, 13 years ago

Got rid of those annoying double_scalar warnings in the windows code (just divide by zero warnings)
File size: 8.3 KB
Line 
1import unittest
2import numpy as num
3
4from urs import calculate_boundary_points, save_boundary_as_urs
5from anuga.coordinate_transforms.redfearn import degminsec2decimal_degrees
6from anuga.geospatial_data.geospatial_data import Geospatial_data
7
8class Test_Urs(unittest.TestCase):
9    def setUp(self):
10        self.verbose = True
11        pass
12
13    def tearDown(self):
14        pass
15       
16    def test_URS_points_needed(self):
17       
18        ll_lat = -21.5
19        ll_long = 114.5
20        grid_spacing = 1./60.
21        lat_amount = 30
22        long_amount = 30
23        zone = 50
24
25        boundary_polygon = [[250000,7660000],[280000,7660000],
26                             [280000,7630000],[250000,7630000]]
27        geo=calculate_boundary_points(boundary_polygon, zone, 
28                              ll_lat, ll_long, grid_spacing, 
29                              lat_amount, long_amount,
30                              verbose=self.verbose)
31        # to test this geo, can info from it be transfered onto the boundary
32        # poly?
33        #Maybe see if I can fit the data to the polygon - have to represent
34        # the poly as points though.
35        #geo.export_points_file("results.txt", as_lat_long=True)
36        results = frozenset(geo.get_data_points(as_lat_long=True))
37        #print 'results',results
38
39        # These are a set of points that have to be in results
40        points = []
41        for i in range(18):
42            lat = -21.0 - 8./60 - grid_spacing * i
43            points.append((lat,degminsec2decimal_degrees(114,35,0))) 
44            points.append((lat,degminsec2decimal_degrees(114,36,0))) 
45            points.append((lat,degminsec2decimal_degrees(114,52,0))) 
46            points.append((lat,degminsec2decimal_degrees(114,53,0)))
47        geo_answer = Geospatial_data(data_points=points,
48                                     points_are_lats_longs=True) 
49        #geo_answer.export_points_file("answer.txt", as_lat_long=True) 
50        answer = frozenset(points)
51       
52        outs = answer.difference(results)
53        #print "outs", outs
54        # This doesn't work.  Though visualising the results shows that
55        # it is correct.
56        #assert answer.issubset(results)
57        # this is why;
58        #point (-21.133333333333333, 114.58333333333333)
59        #result (-21.133333332232368, 114.58333333300342)
60       
61        for point in points:
62            found = False
63            for result in results:
64                if num.allclose(point, result):
65                    found = True
66                    break
67            if not found:
68                assert False
69       
70   
71    def dave_test_URS_points_needed(self):
72        ll_lat = -21.51667
73        ll_long = 114.51667
74        grid_spacing = 2./60.
75        lat_amount = 15
76        long_amount = 15
77
78       
79        boundary_polygon = [[250000,7660000],[280000,7660000],
80                             [280000,7630000],[250000,7630000]]
81        calculate_boundary_points('a_test_example',boundary_polygon,
82                                  ll_lat, ll_long, grid_spacing, 
83                                  lat_amount, long_amount,
84                                  verbose=self.verbose)
85       
86    def X_test_URS_points_neededII(self):
87        ll_lat = -21.5
88        ll_long = 114.5
89        grid_spacing = 1./60.
90        lat_amount = 30
91        long_amount = 30
92
93        # change this so lats and longs are inputed, then converted
94       
95        #boundary_polygon = [[7660000,250000],[7660000,280000],
96        #                     [7630000,280000],[7630000,250000]]
97        calculate_boundary_points(boundary_polygon, ll_lat, ll_long,
98                          grid_spacing, 
99                          lat_amount, long_amount,
100                          verbose=self.verbose)
101       
102    def test_URS_points_northern_hemisphere(self):
103
104        LL_LAT = 8.0
105        LL_LONG = 97.0
106        GRID_SPACING = 2.0/60.0
107        LAT_AMOUNT = 2
108        LONG_AMOUNT = 2
109        ZONE = 47
110
111        #
112        points = []
113        for i in range(2):
114            for j in range(2):
115                points.append((degminsec2decimal_degrees(8,1+i*2,0),
116                               degminsec2decimal_degrees(97,1+i*2,0)))
117        #print "points", points
118        geo_poly = Geospatial_data(data_points=points,
119                                     points_are_lats_longs=True)
120        poly_lat_long = geo_poly.get_data_points(as_lat_long=False,
121                                       isSouthHemisphere=False)
122        #print "seg_lat_long",  poly_lat_long
123       
124      #   geo=URS_points_needed_to_file('test_example_poly3', poly_lat_long,
125#                                   ZONE,
126#                                   LL_LAT, LL_LONG,
127#                                   GRID_SPACING,
128#                                   LAT_AMOUNT, LONG_AMOUNT,
129#                                   isSouthernHemisphere=False,
130#                                   export_csv=True,
131#                                   verbose=self.verbose)
132
133        geo=calculate_boundary_points(poly_lat_long,
134                                  ZONE,
135                                  LL_LAT, LL_LONG,
136                                  GRID_SPACING,
137                                  LAT_AMOUNT, LONG_AMOUNT,
138                                  isSouthHemisphere=False,
139                                  verbose=self.verbose)
140
141        results = frozenset(geo.get_data_points(as_lat_long=True,
142                                                isSouthHemisphere=False))
143
144        #print 'results',results
145
146        # These are a set of points that have to be in results
147        points = [] 
148        for i in range(2):
149            for j in range(2):
150                points.append((degminsec2decimal_degrees(8,i*2,0),
151                               degminsec2decimal_degrees(97,i*2,0)))
152        #print "answer points", points
153        answer = frozenset(points)
154       
155        for point in points:
156            found = False
157            for result in results:
158                if num.allclose(point, result):
159                    found = True
160                    break
161            if not found:
162                assert False
163
164    def covered_in_other_tests_test_URS_points_needed_poly1(self):
165        # Values used for FESA 2007 results
166        # domain in southern hemisphere zone 51       
167        LL_LAT = -50.0
168        LL_LONG = 80.0
169        GRID_SPACING = 2.0/60.0
170        LAT_AMOUNT = 4800
171        LONG_AMOUNT = 3600
172        ZONE = 51
173       
174        poly1 = [[296361.89, 8091928.62],
175                 [429495.07,8028278.82],
176                 [447230.56,8000674.05],
177                 [429661.2,7982177.6],
178                 [236945.9,7897453.16],
179                 [183493.44,7942782.27],
180                 [226583.04,8008058.96]]
181
182        save_boundary_as_urs('test_example_poly2', poly1,
183                                  ZONE,
184                                  LL_LAT, LL_LONG,
185                                  GRID_SPACING,
186                                  LAT_AMOUNT, LONG_AMOUNT,
187                                  verbose=self.verbose)
188       
189
190
191    def covered_in_other_tests_test_URS_points_needed_poly2(self):
192        # Values used for 2004 validation work
193        # domain in northern hemisphere zone 47       
194        LL_LAT = 0.0
195        LL_LONG = 90.0
196        GRID_SPACING = 2.0/60.0
197        LAT_AMOUNT = (15-LL_LAT)/GRID_SPACING
198        LONG_AMOUNT = (100-LL_LONG)/GRID_SPACING
199        ZONE = 47
200       
201        poly2 = [[419336.424,810100.845],
202                 [342405.0281,711455.8026],
203                 [274649.9152,723352.9603],
204                 [272089.092,972393.0131],
205                 [347633.3754,968551.7784],
206                 [427979.2022,885965.2313],
207                 [427659.0993,875721.9386],
208                 [429259.6138,861317.3083],
209                 [436301.8775,840830.723]]
210       
211        save_boundary_as_urs('test_example_poly2', poly2,
212                                  ZONE,
213                                  LL_LAT, LL_LONG,
214                                  GRID_SPACING,
215                                  LAT_AMOUNT, LONG_AMOUNT,
216                                  isSouthernHemisphere=False,
217                                  verbose=self.verbose) 
218       
219
220
221       
222################################################################################
223
224if __name__ == "__main__":
225    suite = unittest.makeSuite(Test_Urs,'test')
226    runner = unittest.TextTestRunner() #verbosity=2)
227    runner.run(suite)
Note: See TracBrowser for help on using the repository browser.