Changeset 3738

Timestamp:

Oct 11, 2006, 10:22:45 AM (18 years ago)

Author:

ole

Message:

Added attributes into clip functions and tested

Location:

anuga_core/source/anuga/geospatial_data

Files:

• geospatial_data.py (modified) (7 diffs)
• test_geospatial_data.py (modified) (13 diffs)

anuga_core/source/anuga/geospatial_data/geospatial_data.py

@@ -5 +5 @@
 
 from os import access, F_OK, R_OK
+from types import DictType
+
 from Numeric import concatenate, array, Float, shape, reshape, ravel, take
 
@@ -124 +126 @@
 
         else:
-###                data = {}    
-                # watch for case where file name and points, attributes etc are provided!!
-                # if file name then all provided info will be removed!
-                self.import_points_file(file_name, delimiter, verbose)
+            # watch for case where file name and points, attributes etc are provided!!
+            # if file name then all provided info will be removed!
+            self.import_points_file(file_name, delimiter, verbose)
                 
-                self.check_data_points(self.data_points)
-                self.set_attributes(self.attributes)
-                self.set_geo_reference(self.geo_reference)
-                self.set_default_attribute_name(default_attribute_name)
+            self.check_data_points(self.data_points)
+            self.set_attributes(self.attributes) 
+            self.set_geo_reference(self.geo_reference)
+            self.set_default_attribute_name(default_attribute_name)
+
+
+        assert self.attributes is None or isinstance(self.attributes, DictType)
+        
 
     def __len__(self):
@@ -159 +164 @@
         """
         
-        from types import DictType
-       
         if attributes is None:
             self.attributes = None
@@ -235 +238 @@
         points = self.get_data_points()    
         inside_indices = inside_polygon(points, polygon, closed)
-        return Geospatial_data(take(points, inside_indices))
+
+        clipped_points = take(points, inside_indices)
+
+        # Clip all attributes
+        attributes = self.get_all_attributes()
+
+        clipped_attributes = {}
+        if attributes is not None:
+            for key, att in attributes.items():
+                clipped_attributes[key] = take(att, inside_indices)
+
+        return Geospatial_data(clipped_points, clipped_attributes)
         
         
@@ -261 +275 @@
         points = self.get_data_points()    
         outside_indices = outside_polygon(points, polygon, closed)
-        return Geospatial_data(take(points, outside_indices))
+
+        clipped_points = take(points, outside_indices)
+
+        # Clip all attributes
+        attributes = self.get_all_attributes()
+
+        clipped_attributes = {}
+        if attributes is not None:
+            for key, att in attributes.items():
+                clipped_attributes[key] = take(att, outside_indices)
+
+        return Geospatial_data(clipped_points, clipped_attributes)
+
     
     def _set_using_lat_long(self,
@@ -350 +376 @@
         """
         Return values for all attributes.
+        The return value is either None or a dictionary (possibly empty).
         """
 
@@ -445 +472 @@
         Note: will throw an IOError if it can't load the file.
         Catch these!
+
+        Post condition: self.attributes dictionary has been set
         """
         

anuga_core/source/anuga/geospatial_data/test_geospatial_data.py

@@ -395 +395 @@
                         [[0.5, 0.5], [1, -0.5], [1.5, 0]])
 
-    def no_test_clip1(self):
+    def test_clip0_with_attributes(self):
+        """test_clip0_with_attributes(self):
+        
+        Test that point sets with attributes can be clipped by a polygon
+        """
+        
+        from anuga.coordinate_transforms.geo_reference import Geo_reference
+        
+        points = [[-1, 4], [0.2, 0.5], [1.0, 2.1], [0.4, 0.3], [3.0, 5.3],
+                  [0, 0], [2.4, 3.3]]
+
+        attributes = [2, -4, 5, 76, -2, 0.1, 3]
+        att_dict = {'att1': attributes,
+                    'att2': array(attributes)+1}
+        
+        G = Geospatial_data(points, att_dict)
+
+        # First try the unit square    
+        U = [[0,0], [1,0], [1,1], [0,1]] 
+        assert allclose(G.clip(U).get_data_points(), [[0.2, 0.5], [0.4, 0.3], [0, 0]])
+        assert allclose(G.clip(U).get_attributes('att1'), [-4, 76, 0.1])
+        assert allclose(G.clip(U).get_attributes('att2'), [-3, 77, 1.1])                
+
+        # Then a more complex polygon
+        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
+        points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
+
+        # This time just one attribute
+        attributes = [2, -4, 5, 76, -2, 0.1]
+        G = Geospatial_data(points, attributes)
+
+        assert allclose(G.clip(polygon).get_data_points(),
+                        [[0.5, 0.5], [1, -0.5], [1.5, 0]])
+        assert allclose(G.clip(polygon).get_attributes(), [-4, 5, 76])
+        
+
+    def test_clip1(self):
         """test_clip1(self):
         
@@ -408 +444 @@
         attributes = [2, -4, 5, 76, -2, 0.1, 3]
         att_dict = {'att1': attributes,
-                    'att2': array(attributes) +1}
+                    'att2': array(attributes)+1}
         G = Geospatial_data(points, att_dict)
         
@@ -417 +453 @@
 
         assert allclose(G.clip(U).get_attributes('att1'), [-4, 76, 0.1])
-        #assert allclose(G.clip(U).get_attributes('att2'), ???)
+        assert allclose(G.clip(U).get_attributes('att2'), [-3, 77, 1.1])                        
         
         # Then a more complex polygon
         points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
-        G = Geospatial_data(points)
+        attributes = [2, -4, 5, 76, -2, 0.1]        
+        G = Geospatial_data(points, attributes)
         polygon = Geospatial_data([[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]])
         
@@ -427 +464 @@
         assert allclose(G.clip(polygon).get_data_points(),
                         [[0.5, 0.5], [1, -0.5], [1.5, 0]])
+        assert allclose(G.clip(polygon).get_attributes(), [-4, 5, 76])
+        
 
     def test_clip0_outside(self):
@@ -438 +477 @@
         points = [[-1, 4], [0.2, 0.5], [1.0, 2.1], [0.4, 0.3], [3.0, 5.3],
                   [0, 0], [2.4, 3.3]]
-        G = Geospatial_data(points)
+        attributes = [2, -4, 5, 76, -2, 0.1, 3]        
+        G = Geospatial_data(points, attributes)
 
         # First try the unit square    
@@ -444 +484 @@
         assert allclose(G.clip_outside(U).get_data_points(),
                         [[-1, 4], [1.0, 2.1], [3.0, 5.3], [2.4, 3.3]])
+        #print G.clip_outside(U).get_attributes()
+        assert allclose(G.clip_outside(U).get_attributes(), [2, 5, -2, 3])        
+        
 
         # Then a more complex polygon
         polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
         points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
-        G = Geospatial_data(points)
+        attributes = [2, -4, 5, 76, -2, 0.1]        
+        G = Geospatial_data(points, attributes)
 
         assert allclose(G.clip_outside(polygon).get_data_points(),
-                        [[0.5, 1.4], [0.5, 1.5], [0.5, -0.5]])    
+                        [[0.5, 1.4], [0.5, 1.5], [0.5, -0.5]])
+        assert allclose(G.clip_outside(polygon).get_attributes(), [2, -2, 0.1])                
 
 
@@ -465 +510 @@
         points = [[-1, 4], [0.2, 0.5], [1.0, 2.1], [0.4, 0.3], [3.0, 5.3],
                   [0, 0], [2.4, 3.3]]
-        G = Geospatial_data(points)
+        attributes = [2, -4, 5, 76, -2, 0.1, 3]        
+        G = Geospatial_data(points, attributes)        
 
         # First try the unit square    
@@ -471 +517 @@
         assert allclose(G.clip_outside(U).get_data_points(),
                         [[-1, 4], [1.0, 2.1], [3.0, 5.3], [2.4, 3.3]])
+        assert allclose(G.clip(U).get_attributes(), [-4, 76, 0.1])        
 
         # Then a more complex polygon
         points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
-        G = Geospatial_data(points)
+        attributes = [2, -4, 5, 76, -2, 0.1]        
+        G = Geospatial_data(points, attributes)
+
         polygon = Geospatial_data([[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]])
         
@@ -480 +529 @@
         assert allclose(G.clip_outside(polygon).get_data_points(),
                         [[0.5, 1.4], [0.5, 1.5], [0.5, -0.5]])
+        assert allclose(G.clip_outside(polygon).get_attributes(), [2, -2, 0.1])
+        
 
 
@@ -493 +544 @@
         points = [[-1, 4], [0.2, 0.5], [1.0, 2.1], [0.4, 0.3], [3.0, 5.3],
                   [0, 0], [2.4, 3.3]]
-
-        G = Geospatial_data(points)
+        attributes = [2, -4, 5, 76, -2, 0.1, 3]        
+        G = Geospatial_data(points, attributes)
 
         # First try the unit square    
@@ -500 +551 @@
         G1 = G.clip(U)
         assert allclose(G1.get_data_points(),[[0.2, 0.5], [0.4, 0.3], [0, 0]])
+        assert allclose(G.clip(U).get_attributes(), [-4, 76, 0.1])
         
         G2 = G.clip_outside(U)
         assert allclose(G2.get_data_points(),[[-1, 4], [1.0, 2.1],
                                               [3.0, 5.3], [2.4, 3.3]])
+        assert allclose(G.clip_outside(U).get_attributes(), [2, 5, -2, 3])                
 
         
@@ -510 +563 @@
                      [[-1, 4], [1.0, 2.1], [3.0, 5.3], [2.4, 3.3]]
 
+        new_attributes = [-4, 76, 0.1, 2, 5, -2, 3]                 
+        
         assert allclose((G1+G2).get_data_points(), new_points)
+        assert allclose((G1+G2).get_attributes(), new_attributes)
 
         G = G1+G2
@@ -523 +579 @@
         # Check result
         assert allclose(G3.get_data_points(), new_points)        
-        
+        assert allclose(G3.get_attributes(), new_attributes)        
         
         os.remove(FN)