Changeset 4223

Timestamp:

Feb 6, 2007, 2:28:47 PM (18 years ago)

Author:

duncan

Message:

Start of improved urs2sww, along with other stuff.

Location:

anuga_core/source/anuga

Files:

• coordinate_transforms/redfearn.py (modified) (1 diff)
• fit_interpolate/benchmark_least_squares.py (modified) (2 diffs)
• fit_interpolate/run_long_benchmark.py (modified) (3 diffs)
• geospatial_data/geospatial_data.py (modified) (3 diffs)
• geospatial_data/test_geospatial_data.py (modified) (6 diffs)
• shallow_water/data_manager.py (modified) (2 diffs)
• shallow_water/test_data_manager.py (modified) (1 diff)

anuga_core/source/anuga/coordinate_transforms/redfearn.py

@@ -189 +189 @@
 
     Assume the false_easting and false_northing are the same for each list.
-    If points end up belonging to different UTM zones, an ANUGAerror is thrown.    
+    If points end up in different UTM zones, an ANUGAerror is thrown.    
     """
 

anuga_core/source/anuga/fit_interpolate/benchmark_least_squares.py

@@ -68 +68 @@
               use_least_squares=False,
               use_file_type=None,
-              save=False):
+              save=False,
+              verbose=False):
         '''
         num_of_points 
@@ -117 +118 @@
                                          points_dict['point_attributes'])
                     G1.export_points_file(fileName, absolute=True)
-                    calc = interp.fit(fileName)
+                    calc = interp.fit(fileName), verbose=verbose)
                     os.remove(fileName)
                     

anuga_core/source/anuga/fit_interpolate/run_long_benchmark.py

@@ -19 +19 @@
 #num_of_points_list = [10]
 #maxArea_list = [0.1, 0.001]
-num_of_points_list = [10, 100] #, 10000, 100000] #, 10000000]
-maxArea_list = [0.1, 0.001] #, 0.00001, 0.0000001]
+num_of_points_list = [10000, 100000] #, 10000, 100000] #, 10000000]
+maxArea_list = [0.1, 0.01, 0.001] #, 0.00001, 0.0000001]
 max_points_per_cell_list = [8]
-use_file_type_list = [None,'txt']
+use_file_type_list = [None,'txt', 'pts']
 
 fd = open(ofile,'a')
@@ -37 +37 @@
 
 
-for use_file_type in use_file_type_list:
-    for is_fit in is_fit_list:
-        for num_of_points in num_of_points_list:
-            for maxArea in maxArea_list:
+for maxArea in maxArea_list:
+    for use_file_type in use_file_type_list:
+        for is_fit in is_fit_list:
+            for num_of_points in num_of_points_list:
                 for max_points_per_cell in max_points_per_cell_list:
     
@@ -52 +52 @@
                     print "mem", mem
                     fd.write(str(use_file_type) + delimiter +
-                             str(is_fit) + delimiter +
                              str(num_of_points) + delimiter +
                              str(maxArea) + delimiter +
                              str(num_tri) + delimiter +
                              str(max_points_per_cell) + delimiter +
+                             str(is_fit) + delimiter +
                              str(mem)  + delimiter +
                              str(time) + delimiter + "\n")

anuga_core/source/anuga/geospatial_data/geospatial_data.py

@@ -14 +14 @@
 
 from Scientific.IO.NetCDF import NetCDFFile    
-    
+from anuga.coordinate_transforms.lat_long_UTM_conversion import UTMtoLL    
 from anuga.utilities.numerical_tools import ensure_numeric
-from anuga.coordinate_transforms.geo_reference import Geo_reference, TitleError
+from anuga.coordinate_transforms.geo_reference import Geo_reference, \
+     TitleError, DEFAULT_ZONE
 from anuga.coordinate_transforms.redfearn import convert_from_latlon_to_utm
 from anuga.utilities.anuga_exceptions import ANUGAError
@@ -332 +333 @@
         return self.geo_reference
        
-    def get_data_points(self, absolute=True, geo_reference=None):
+    def get_data_points(self, absolute=True, geo_reference=None,
+                        as_lat_long=False):
         """Get coordinates for all data points as an Nx2 array
 
@@ -344 +346 @@
         Default: absolute is True.
         """
-
+        if as_lat_long is True:
+            msg = "Points need a zone to be converted into lats and longs"
+            assert self.geo_reference is not None, msg
+            zone = self.geo_reference.get_zone()
+            assert self.geo_reference.get_zone() is not DEFAULT_ZONE, msg
+            lats_longs = []
+            for point in self.get_data_points(True):
+                ### UTMtoLL(northing, easting, zone,
+                lat_calced, long_calced = UTMtoLL(point[1],point[0], zone)
+                lats_longs.append([lat_calced, long_calced])
+            return lats_longs
+            
         if absolute is True and geo_reference is None:
             return self.geo_reference.get_absolute(self.data_points)

anuga_core/source/anuga/geospatial_data/test_geospatial_data.py

@@ -7 +7 @@
 from math import sqrt, pi
 import tempfile
+from sets import ImmutableSet
 
 from anuga.geospatial_data.geospatial_data import *
@@ -162 +163 @@
         assert allclose(results, points_rel)
 
-        
+  
+    def test_get_data_points_lat_long(self):
+        # lat long [-30.],[130]
+        #Zone:   52    
+        #Easting:  596450.153  Northing: 6680793.777 
+        # lat long [-32.],[131]
+        #Zone:   52    
+        #Easting:  688927.638  Northing: 6457816.509 
+        
+        points_Lat_long = [[-30.,130], [-32,131]]
+        
+        spatial = Geospatial_data(latitudes=[-30, -32.],
+                                  longitudes=[130, 131])
+
+        results = spatial.get_data_points(as_lat_long=True)
+        #print "test_get_data_points_lat_long - results", results
+        #print "points_Lat_long",points_Lat_long 
+        assert allclose(results, points_Lat_long)
+      
+    def test_get_data_points_lat_longII(self):
+        # x,y  North,east long,lat
+        boundary_polygon = [[ 250000, 7630000]]
+        zone = 50
+        
+        geo_reference = Geo_reference(zone=zone)
+        geo = Geospatial_data(boundary_polygon,geo_reference=geo_reference)
+        seg_lat_long = geo.get_data_points(as_lat_long=True)
+        lat_result = degminsec2decimal_degrees(-21,24,54)
+        long_result = degminsec2decimal_degrees(114,35,17.89)
+        #print "seg_lat_long", seg_lat_long [0][0]
+        #print "lat_result",lat_result 
+        assert allclose(seg_lat_long[0][0], lat_result)#lat
+        assert allclose(seg_lat_long[0][1], long_result)#long
+
+              
     def test_set_geo_reference(self):
         points_ab = [[12.5,34.7],[-4.5,-60.0]]
@@ -1058 +1093 @@
 
     def verbose_test_load_pts_blocking(self):
-        """ test_load_csv(self):
-        space delimited
-        """
+        
         import os
        
@@ -1130 +1163 @@
                         [10.0, 0.0])
         assert allclose(geo_list[1].get_data_points(),
-                        [[1.0, 0.0]])        
+                        [[1.0, 0.0],[0.0, 1.0] ])        
         assert allclose(geo_list[1].get_attributes(attribute_name='elevation'),
-                        [10.4])
+                        [10.0, 0.0])
            
         os.remove(fileName)  
@@ -1981 +2014 @@
             self.failUnless(0 ==1,  'Error not thrown error!')
 
+
+    def test_lat_long_set(self):
+        lat_gong = degminsec2decimal_degrees(-34,30,0.)
+        lon_gong = degminsec2decimal_degrees(150,55,0.)
+        
+        lat_2 = degminsec2decimal_degrees(-34,00,0.)
+        lon_2 = degminsec2decimal_degrees(150,00,0.)
+        p1 = (lat_gong, lon_gong)
+        p2 = (lat_2, lon_2)
+        points = ImmutableSet([p1, p2, p1])
+        gsd = Geospatial_data(data_points=list(points),
+                              points_are_lats_longs=True)
+
+        points = gsd.get_data_points(absolute=True)
+        
+        assert allclose(points[0][0], 308728.009)
+        assert allclose(points[0][1], 6180432.601)
+        assert allclose(points[1][0],  222908.705)
+        assert allclose(points[1][1], 6233785.284)
+        self.failUnless(gsd.get_geo_reference().get_zone() == 56,
+                        'Bad zone error!')
+        points = gsd.get_data_points(as_lat_long=True)
+        #print "test_lat_long_set points", points
+        assert allclose(points[1][0], -34)
+        assert allclose(points[1][1], 150)
+
     def test_len(self):
         
@@ -2025 +2084 @@
 if __name__ == "__main__":
 
-#    suite = unittest.makeSuite(Test_Geospatial_data, 'test_split')
-#    suite = unittest.makeSuite(Test_Geospatial_data, 'test_clip0')
-    suite = unittest.makeSuite(Test_Geospatial_data, 'test')
+    suite = unittest.makeSuite(Test_Geospatial_data, 'test_lat_long_set')
+    #suite = unittest.makeSuite(Test_Geospatial_data, 'verbose_test_load_pts_blocking')
+    #suite = unittest.makeSuite(Test_Geospatial_data, 'test')
     runner = unittest.TextTestRunner()
     runner.run(suite)

anuga_core/source/anuga/shallow_water/data_manager.py

@@ -68 +68 @@
 
 
+from anuga.coordinate_transforms.redfearn import redfearn
 from anuga.coordinate_transforms.geo_reference import Geo_reference
 from anuga.geospatial_data.geospatial_data import Geospatial_data
@@ -4404 +4405 @@
     return long, lat, quantity
 
-    ####  END URS 2 SWW  ###     
+    ####  END URS 2 SWW  ###
+
+    #### URS UNGRIDDED 2 SWW ###
+
+def URS_points_needed_to_file(boundary_polygon, ll_lat, ll_long,
+                              grid_spacing, 
+                      lat_amount, long_amount, zone=None):
+    
+    geo = URS_points_needed(boundary_polygon, ll_lat, ll_long, grid_spacing, 
+                      lat_amount, long_amount, zone=zone)
+    
+def URS_points_needed(boundary_polygon, ll_lat, ll_long, grid_spacing, 
+                      lat_amount, long_amount, zone=None):
+    """
+
+    boundary_polygon - a list of points that describes a polygon.
+                      The last point is assumed ot join the first point.
+                      This is in UTM (lat long would be better though)
+
+    ll_lat - lower left latitude, in decimal degrees
+    ll-long - lower left longitude, in decimal degrees
+    grid_spacing - in deciamal degrees
+
+    """
+    from sets import ImmutableSet
+    
+    msg = "grid_spacing can not be zero"
+    assert not grid_spacing ==0, msg 
+    a = boundary_polygon
+    # List of segments.  Each segment is two points.
+    segs = [i and [a[i-1], a[i]] or [a[len(a)-1], a[0]] for i in range(len(a))]
+
+    # convert the segs to Lat's and longs.
+    if zone is None:
+        # Assume the zone of the segments is the same as the lower left
+        # corner of the lat long data
+        zone, _, _ = redfearn(ll_lat, ll_long)
+    lat_long_set = ImmutableSet()
+    for seg in segs:
+        points_lat_long = points_needed(seg, ll_lat, ll_long, grid_spacing, 
+                      lat_amount, long_amount, zone)
+        lat_long_set |= ImmutableSet(points_lat_long)
+    #print "lat_long_set",lat_long_set 
+    geo = Geospatial_data(data_points=list(lat_long_set),
+                              points_are_lats_longs=True)
+    return geo
+    
+def points_needed(seg, ll_lat, ll_long, grid_spacing, 
+                  lat_amount, long_amount, zone):
+    """
+    return a list of the points, in lats and longs that are needed to
+    interpolate any point on the segment.
+    """
+    #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)
+    #print "seg_lat_long", seg_lat_long
+    # 1.415 = 2^0.5, rounded up....
+    buffer = 1.415 * grid_spacing
+    
+    #
+    
+    max_lat = max(seg_lat_long[0][0], seg_lat_long[1][0]) + buffer
+    max_long = max(seg_lat_long[0][1], seg_lat_long[1][1]) + buffer
+    min_lat = max(seg_lat_long[0][0], seg_lat_long[1][0]) - buffer
+    min_long = max(seg_lat_long[0][1], seg_lat_long[1][1]) - buffer
+
+    #print "ll_lat", ll_lat
+    #print "ll_long", ll_long
+    #print "max_lat", max_lat
+    #print "max_long", max_long 
+    #print "min_lat", min_lat
+    #print "min_long", min_long
+    first_row = (min_long - ll_long)/grid_spacing
+    # To round up
+    first_row_long = int(round(first_row + 0.5))
+    #print "first_row", first_row_long
+
+    last_row = (max_long - ll_long)/grid_spacing # round down
+    last_row_long = int(round(last_row))
+    #print "last_row",last_row _long
+
+    
+    first_row = (min_lat - ll_lat)/grid_spacing
+    # To round up
+    first_row_lat = int(round(first_row + 0.5))
+    #print "first_row", first_row_lat
+
+    last_row = (max_lat - ll_lat)/grid_spacing # round down
+    last_row_lat = int(round(last_row))
+    #print "last_row",last_row_lat
+
+    points_lat_long = []
+    # Create a list of the lat long points to include.
+    for index_lat in range(first_row_lat, last_row_lat + 1):
+        for index_long in range(first_row_long, last_row_long + 1):
+            lat = ll_lat + index_lat*grid_spacing
+            long = ll_long + index_long*grid_spacing
+            points_lat_long.append((lat, long)) #must be hashable
+    #print "points_lat_long", points_lat_long
+    return points_lat_long    
+    
+    #### END URS UNGRIDDED 2 SWW ###
+
 #-------------------------------------------------------------
 if __name__ == "__main__":

anuga_core/source/anuga/shallow_water/test_data_manager.py

@@ -5035 +5035 @@
             raise msg
         
-    def trial_loading(self):
-        basename_in = 'karratha'
-        basename_out = basename_in
-        urs2sww(basename_in, basename_out, remove_nc_files=True,
-                zscale=10000000)
     #### END TESTS FOR URS 2 SWW  ###
 
+    #### TESTS URS UNGRIDDED 2 SWW ###
+    def test_URS_points_needed(self):
+        ll_lat = -21.5
+        ll_long = 114.5
+        grid_spacing = 1./60.
+        lat_amount = 30
+        long_amount = 30
+
+        boundary_polygon = [[250000,7660000],[280000,7660000],
+                             [280000,7630000],[250000,7630000]]
+        geo=URS_points_needed(boundary_polygon,
+                              ll_lat, ll_long, grid_spacing, 
+                              lat_amount, long_amount)
+        # to test this geo, can info from it be transfered onto the boundary
+        # poly?
+        #Maybe see if I can fit the data to the polygon - have to represent
+        # the poly as points though.
+        
+        results = geo.get_data_points(as_lat_long=True)
+        #print 'results',results
+        
+    def X_test_URS_points_neededII(self):
+        ll_lat = -21.5
+        ll_long = 114.5
+        grid_spacing = 1./60.
+        lat_amount = 30
+        long_amount = 30
+
+        # change this so lats and longs are inputed, then converted
+        #boundary_polygon = [[7660000,250000],[7660000,280000],
+        #                     [7630000,280000],[7630000,250000]]
+        URS_points_needed(boundary_polygon, ll_lat, ll_long, grid_spacing, \
+                      lat_amount, long_amount)
+        
+    #### END TESTS URS UNGRIDDED 2 SWW ###
         
 #-------------------------------------------------------------
 if __name__ == "__main__":
-    #suite = unittest.makeSuite(Test_Data_Manager,'test_urs2sww_m')
+    #suite = unittest.makeSuite(Test_Data_Manager,'test_URS_points_needed')
     #suite = unittest.makeSuite(Test_Data_Manager,'test_get_min_max_indexes_lat_ascending')
     #suite = unittest.makeSuite(Test_Data_Manager,'test_ferret2sww_lat_long')