Changeset 2590

Timestamp:

Mar 24, 2006, 2:43:10 PM (19 years ago)

Author:

ole

Message:

Clean up

Location:

inundation/geospatial_data

Files:

• geospatial_data.py (modified) (4 diffs)
• test_geospatial_data.py (modified) (2 diffs)

inundation/geospatial_data/geospatial_data.py

@@ -169 +169 @@
        
     def get_data_points(self, absolute = True):
+        """Get coordinates for all data points as an Nx2 array
+
+        If absolute is True absolute UTM coordinates are returned otherwise
+        returned coordinates are relative to the internal georeference's
+        xll and yll corners.
+
+        Default: absolute is True.
+        """
         if absolute is True:
             xll = self.geo_reference.xllcorner
             yll = self.geo_reference.yllcorner
             return self.data_points + [xll, yll]
-            
         else: 
             return self.data_points
@@ -224 +231 @@
         
         if zone1 == zone2:
-            a_rel_points = self.get_data_points()
-            b_rel_points = other.get_data_points()
-            
-            # subtracts xll and yll from self's and other's
-            # reletive data point and concatenates
-            c_points = concatenate((a_rel_points-[xll, yll],
-                                    b_rel_points-[xll, yll]), axis = 0)
+            a_absolute_points = self.get_data_points(absolute=True)
+            b_absolute_points = other.get_data_points(absolute=True)
+            
+            # subtract xll and yll from self's and other's
+            # absolute data points and concatenate
+            c_points = concatenate((a_absolute_points-[xll, yll],
+                                    b_absolute_points-[xll, yll]),
+                                   axis = 0)
 
             new_attributes = {}
@@ -388 +396 @@
     line = fd.readline()
     numbers = clean_line(line,delimiter)
-    print 'read xya numbers', numbers
+    #print 'read xya numbers', numbers
     while len(numbers) > 1:
         if numbers != []:
@@ -588 +596 @@
     that are shared. This will require some work and maybe __subtract__ function 
     '''
-#    print 'load %s' %add_file1
-#    if access(add_file1,F_OK) == 1 :
-#        print'file1', add_file1
-
-#    print'file1 again', add_file1
-#    print'file2 again', add_file2
-    print 'create G1'
     G1 = Geospatial_data(file_name = add_file1)
-    
-#    print 'G1 dict'
- #   G1_points_dict = geospatial_data2points_dictionary(G1)
-
-#    print 'G1 xya file'
-#    export_points_file("v:\inundation\onslow_tsunami_scenario_2006\\topographies\onshore.xya", G1_points_dict)
-
-    print 'create G2'
-#    print 'load %s' %add_file2
-#    if access(add_file2,F_OK) == 1 :
-#        print'file2', add_file2
     G2 = Geospatial_data(file_name = add_file2)
-    print 'G2 dict'
-#    G2_points_dict = geospatial_data2points_dictionary(G2)
-
     new_add_file2 = add_file2[:-4] + '.pts' 
-    print 'G2 pts file', new_add_file2
-    
-#    export_points_file(new_add_file2, G2_points_dict)
-    print 'G1+g2'    
+
     G = G1 + G2
+    
     #FIXME remove dependance on points to dict in export only!
     G_points_dict = geospatial_data2points_dictionary(G)
-    
     export_points_file(results_file, G_points_dict)
     

inundation/geospatial_data/test_geospatial_data.py

@@ -202 +202 @@
                                               [1.0, 2.1], [3.0, 5.3]])
         
-    def test_add1 (self):
+    def test_add_with_geo (self):
         """
         Difference in Geo_reference resolved
@@ -208 +208 @@
         points1 = [[1.0, 2.1], [3.0, 5.3]]
         points2 = [[5.0, 6.1], [6.0, 3.3]]
-        attributes1= [2, 4]
+        attributes1 = [2, 4]
         attributes2 = [5, 76]
         geo_ref1= Geo_reference(55, 1.0, 2.0)
-        geo_ref2 = Geo_reference(zone = 55, xllcorner = 0.1,
-                                 yllcorner = 3.0, datum = 'wgs84',
-                                 projection = 'UTM', units = 'm')
+        geo_ref2 = Geo_reference(zone=55,
+                                 xllcorner=0.1,
+                                 yllcorner=3.0,
+                                 datum='wgs84',
+                                 projection='UTM',
+                                 units='m')
                                 
         G1 = Geospatial_data(points1, attributes1, geo_ref1)
         G2 = Geospatial_data(points2, attributes2, geo_ref2)
 
+        #Check that absolute values are as expected
+        P1 = G1.get_data_points(absolute=True)
+        assert allclose(P1, [[2.0, 4.1], [4.0, 7.3]])
+
+        P2 = G2.get_data_points(absolute=True)
+        assert allclose(P2, [[5.1, 9.1], [6.1, 6.3]])        
+        
         G = G1 + G2
         
         assert allclose(G.get_geo_reference().get_xllcorner(), 0.1)
         assert allclose(G.get_geo_reference().get_yllcorner(), 2.0)
-#        print 'G data points=', G.get_data_points()
-        assert allclose(G.get_data_points(), [[2.0, 4.1], [4.0, 7.3], [5.1, 9.1], [6.1, 6.3]])                             
-
+
+        P = G.get_data_points(absolute=True)
+        assert allclose(P, concatenate( (P1,P2) ))
+
+        P_relative = G.get_data_points(absolute=False)
+        assert allclose(P_relative, P - [0.1, 2.0])                      
+        assert allclose(P, [[2.0, 4.1], [4.0, 7.3],
+                            [5.1, 9.1], [6.1, 6.3]])