Changeset 5254

Timestamp:

Apr 30, 2008, 4:34:19 PM (17 years ago)

Author:

duncan

Message:

Adding northern hemisphere tests for URS_points_needed.

Location:

anuga_core/source/anuga/shallow_water

Files:

• data_manager.py (modified) (7 diffs)
• test_data_manager.py (modified) (3 diffs)

anuga_core/source/anuga/shallow_water/data_manager.py

@@ -4439 +4439 @@
     ll-long - lower left longitude, in decimal degrees
     grid_spacing - in deciamal degrees
-    lat_amount - number of latitudes
-    long_amount- number of longs
+    lat_amount - number of latitudes 
+    long_amount- number of longs 
 
 
@@ -4456 +4456 @@
             file_name = file_name[:-4] + ".csv"
         geo.export_points_file(file_name)
+    return geo
 
 def URS_points_needed(boundary_polygon, zone, ll_lat,
@@ -4531 +4532 @@
                   isSouthHemisphere):
     """
-    seg is one point, in UTM
+    seg is two points, in UTM
     return a list of the points, in lats and longs that are needed to
     interpolate any point on the segment.
@@ -4538 +4539 @@
     #print "zone",zone 
     geo_reference = Geo_reference(zone=zone)
+    #print "seg", seg
     geo = Geospatial_data(seg,geo_reference=geo_reference)
     seg_lat_long = geo.get_data_points(as_lat_long=True,
                                        isSouthHemisphere=isSouthHemisphere)
+    #print "seg_lat_long", seg_lat_long
     # 1.415 = 2^0.5, rounded up....
     sqrt_2_rounded_up = 1.415
@@ -4550 +4553 @@
     min_long = min(seg_lat_long[0][1], seg_lat_long[1][1]) - buffer
 
+    #print "min_long", min_long
+    #print "ll_long", ll_long
+    #print "grid_spacing", grid_spacing
     first_row = (min_long - ll_long)/grid_spacing
     # To round up
@@ -4577 +4583 @@
     for index_lat in range(first_row_lat, last_row_lat + 1):
         for index_long in range(first_row_long, last_row_long + 1):
+            #print "index_lat", index_lat
+            #print "index_long", index_long
             lat = ll_lat + index_lat*grid_spacing
             long = ll_long + index_long*grid_spacing
@@ -4598 +4606 @@
     x2 = seg[1][0]
     y2 = seg[1][1]
-
     x2_1 = x2-x1
     y2_1 = y2-y1
-    d = abs((x2_1)*(y1-y0)-(x1-x0)*(y2_1))/sqrt((x2_1)*(x2_1)+(y2_1)*(y2_1))
+    try:
+        d = abs((x2_1)*(y1-y0)-(x1-x0)*(y2_1))/sqrt( \
+            (x2_1)*(x2_1)+(y2_1)*(y2_1))
+    except ZeroDivisionError:
+        #print "seg", seg
+        #print "x0", x0
+        #print "y0", y0
+        if sqrt((x2_1)*(x2_1)+(y2_1)*(y2_1)) == 0 and \
+           abs((x2_1)*(y1-y0)-(x1-x0)*(y2_1)) == 0:
+            return True
+        else:
+            return False
+    
     if d <= max_distance:
         return True

anuga_core/source/anuga/shallow_water/test_data_manager.py

@@ -6117 +6117 @@
                           lat_amount, long_amount,
                           verbose=self.verbose)
-
-    def test_URS_points_needed_poly1(self):
+        
+    def test_URS_points_northern_hemisphere(self):
+               
+        LL_LAT = 8.0
+        LL_LONG = 97.0
+        GRID_SPACING = 2.0/60.0
+        LAT_AMOUNT = 2
+        LONG_AMOUNT = 2
+        ZONE = 47
+
+        # 
+        points = []
+        for i in range(2):
+            for j in range(2):
+                points.append((degminsec2decimal_degrees(8,1+i*2,0),
+                               degminsec2decimal_degrees(97,1+i*2,0)))
+        #print "points", points
+        geo_poly = Geospatial_data(data_points=points,
+                                     points_are_lats_longs=True)
+        poly_lat_long = geo_poly.get_data_points(as_lat_long=False,
+                                       isSouthHemisphere=False)
+        #print "seg_lat_long",  poly_lat_long
+        
+        #geo=URS_points_needed_to_file('test_example_poly3', poly_lat_long,
+        geo=URS_points_needed(poly_lat_long,
+                                  ZONE,
+                                  LL_LAT, LL_LONG,
+                                  GRID_SPACING,
+                                  LAT_AMOUNT, LONG_AMOUNT,
+                                  isSouthHemisphere=False,
+                                  verbose=self.verbose) 
+        results = ImmutableSet(geo.get_data_points(as_lat_long=True,
+                                  isSouthHemisphere=False))
+        #print 'results',results
+
+        # These are a set of points that have to be in results
+        points = [] 
+        for i in range(2):
+            for j in range(2):
+                points.append((degminsec2decimal_degrees(8,i*2,0),
+                               degminsec2decimal_degrees(97,i*2,0)))
+        #print "answer points", points
+        answer = ImmutableSet(points)
+        
+        for point in points:
+            found = False
+            for result in results:
+                if allclose(point, result):
+                    found = True
+                    break
+            if not found:
+                assert False
+        
+
+    def covered_in_other_tests_test_URS_points_needed_poly1(self):
         # Values used for FESA 2007 results
         # domain in southern hemisphere zone 51        
@@ -6141 +6194 @@
                                   GRID_SPACING,
                                   LAT_AMOUNT, LONG_AMOUNT,
-                                  verbose=self.verbose)         
-
-
-    def test_URS_points_needed_poly2(self):
+                                  verbose=self.verbose)
+        
+
+
+    def covered_in_other_tests_test_URS_points_needed_poly2(self):
         # Values used for 2004 validation work
         # domain in northern hemisphere zone 47        
@@ -7562 +7616 @@
     suite = unittest.makeSuite(Test_Data_Manager,'test')
     #suite = unittest.makeSuite(Test_Data_Manager,'test_get_flow_through_cross_section')
-    #suite = unittest.makeSuite(Test_Data_Manager,'Xtest')
+    #suite = unittest.makeSuite(Test_Data_Manager,'covered_')