Changeset 3969


Ignore:
Timestamp:
Nov 13, 2006, 1:53:33 PM (17 years ago)
Author:
nick
Message:

added split and get sample modules to geospatial_data
also removed xxx_add_data_points this was not being used and is old

Location:
anuga_core/source/anuga/geospatial_data
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/geospatial_data/geospatial_data.py

    r3745 r3969  
    77from types import DictType
    88
    9 from Numeric import concatenate, array, Float, shape, reshape, ravel, take
     9from Numeric import concatenate, array, Float, shape, reshape, ravel, take, \
     10                        size, shape
     11from RandomArray import randint
     12#from MA import tolist
    1013
    1114from anuga.utilities.numerical_tools import ensure_numeric
     
    223226          regarded as belonging to the polygon (closed = True)
    224227          or not (closed = False). Default is True.
    225          
     228           
    226229        Output
    227230          New geospatial data object representing points inside
     
    239242        inside_indices = inside_polygon(points, polygon, closed)
    240243
    241         clipped_points = take(points, inside_indices)
     244        clipped_G = self.get_sample(inside_indices)
     245#        clipped_points = take(points, inside_indices)
    242246
    243247        # Clip all attributes
    244         attributes = self.get_all_attributes()
    245 
    246         clipped_attributes = {}
    247         if attributes is not None:
    248             for key, att in attributes.items():
    249                 clipped_attributes[key] = take(att, inside_indices)
    250 
    251         return Geospatial_data(clipped_points, clipped_attributes)
     248#        attributes = self.get_all_attributes()
     249
     250#        clipped_attributes = {}
     251#        if attributes is not None:
     252#            for key, att in attributes.items():
     253#                clipped_attributes[key] = take(att, inside_indices)
     254
     255#        return Geospatial_data(clipped_points, clipped_attributes)
     256        return clipped_G
    252257       
    253258       
     
    276281        outside_indices = outside_polygon(points, polygon, closed)
    277282
    278         clipped_points = take(points, outside_indices)
     283        clipped_G = self.get_sample(outside_indices)
     284
     285#        clipped_points = take(points, outside_indices)
    279286
    280287        # Clip all attributes
    281         attributes = self.get_all_attributes()
    282 
    283         clipped_attributes = {}
    284         if attributes is not None:
    285             for key, att in attributes.items():
    286                 clipped_attributes[key] = take(att, outside_indices)
    287 
    288         return Geospatial_data(clipped_points, clipped_attributes)
     288#        attributes = self.get_all_attributes()
     289
     290#        clipped_attributes = {}
     291#        if attributes is not None:
     292#            for key, att in attributes.items():
     293#                clipped_attributes[key] = take(att, outside_indices)
     294
     295#        return Geospatial_data(clipped_points, clipped_attributes)
     296        return clipped_G
    289297
    290298   
     
    565573            msg = 'Unknown file type %s ' %file_name
    566574            raise IOError, msg
    567    
    568 
     575       
     576    def get_sample(self, indices):
     577        """ Returns a object which is a subset of the original
     578        and the data points and attributes in this new object refer to
     579        the indices provided
     580       
     581        Input
     582            indices- a list of integers that represent the new object
     583        Output
     584            New geospatial data object representing points specified by
     585            the indices
     586            """
     587        #FIXME: add the geo_reference to this
     588       
     589        points = self.get_data_points()
     590        sampled_points = take(points, indices)
     591
     592        attributes = self.get_all_attributes()
     593
     594        sampled_attributes = {}
     595        if attributes is not None:
     596            for key, att in attributes.items():
     597                sampled_attributes[key] = take(att, indices)
     598
     599        return Geospatial_data(sampled_points, sampled_attributes)
     600   
     601   
     602    def split(self, factor=0.5):
     603        """Returns two geospatial_data object, first is size of the 'factor'
     604        smaller the original one and the second is the remainer. The two new
     605        object are disjoin set of each other.
     606       
     607        Points of the two new object have selected RANDOMLY.
     608        AND if factor is a decimal it will round (2.25 to 2 and 2.5 to 3)
     609       
     610       
     611        Input - the factor which to split the object, if 0.1 then 10% of the
     612            object will be returned
     613       
     614        Output - two geospatial_data objects that are disjoint sets of the
     615            original
     616            """
     617       
     618        i=0
     619        self_size = len(self)
     620#        print 'size', self_size
     621        random_list = []
     622        remainder_list = []
     623        new_size = round(factor*self_size)
     624   #     print'Split original %s by %s' %(self_size, factor)
     625   #     print'New samples are %s and %s in size' %(int(round(factor*self_size)),int(self_size-new_size))
     626       
     627        #find unique random numbers
     628        while i < new_size:
     629            random_num = randint(0,self_size)
     630            if random_num not in random_list:
     631                random_list.append(random_num)
     632                i=i+1
     633
     634        #Make list of opposite to random_list
     635        for i in range(0,self_size,1):
     636            remainder_list.append(i)
     637#        print 'start remainer',remainder_list
     638
     639        #need to sort and reverse so the pop() works correctly
     640        random_list.sort()
     641        random_list.reverse()
     642        for i in random_list:
     643            remainder_list.pop(i)
     644
     645        #get new samples
     646        G1 = self.get_sample(random_list)
     647        G2 = self.get_sample(remainder_list)
     648
     649        return G1, G2
    569650
    570651
     
    830911    return numbers
    831912           
    832 def xxx_add_points_files(add_file1, add_file2, results_file):
    833     """ adds the points and attruibutes of 2 pts or xya files and
    834     writes it to a pts file
    835    
    836     NOTE will add the function to check and remove points from one set
    837     that are shared. This will require some work and maybe __subtract__ function
    838     """
    839    
    840     G1 = Geospatial_data(file_name = add_file1)
    841     G2 = Geospatial_data(file_name = add_file2)
    842     new_add_file2 = add_file2[:-4] + '.pts'
    843 
    844     G = G1 + G2
    845    
    846     #FIXME remove dependance on points to dict in export only!
    847 #    G_points_dict = geospatial_data2points_dictionary(G)
    848 #    export_points_file(results_file, G_points_dict)
    849 
    850 #    G_points_dict = geospatial_data2points_dictionary(G)
    851 
    852     G.export_points_file(results_file)
    853 
    854 # '
    855913def ensure_absolute(points, geo_reference = None):
    856914    """
     
    924982    points = Geospatial_data(data_points=points, geo_reference=geo)       
    925983    return points
     984
     985#def file2xya(filename):
     986   
     987#    G = Geospatial_data(filename)
     988#    G.export_points_file()
     989
    926990             
  • anuga_core/source/anuga/geospatial_data/test_geospatial_data.py

    r3797 r3969  
    16271627        G = Geospatial_data(points)
    16281628        self.failUnless(4 ==len(G),  'Len error!')
     1629       
     1630    def test_split(self):
     1631        """test if the results from spilt are disjoin sets"""
     1632       
     1633        points = [[1.0, 1.0], [1.0, 2.0],[1.0, 3.0], [1.0, 4.0], [1.0, 5.0],
     1634                  [2.0, 1.0], [2.0, 2.0],[2.0, 3.0], [2.0, 4.0], [2.0, 5.0],
     1635                  [3.0, 1.0], [3.0, 2.0],[3.0, 3.0], [3.0, 4.0], [3.0, 5.0],
     1636                  [4.0, 1.0], [4.0, 2.0],[4.0, 3.0], [4.0, 4.0], [4.0, 5.0],
     1637                  [5.0, 1.0], [5.0, 2.0],[5.0, 3.0], [5.0, 4.0], [5.0, 5.0]]
     1638        attributes = {'depth':[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,
     1639                      14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25],'speed':[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,
     1640                      14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25]}
     1641        G = Geospatial_data(points, attributes)
     1642#        print G.get_data_points()
     1643#        print G.get_attributes()
     1644
     1645        factor = 0.22
     1646
     1647        G1, G2  = G.split(factor) #will return G1 with 10% for points and G2 with 90%
     1648       
     1649#        print 'len(G): %s  len(G1): %s len(G2): %s' %(len(G), len(G1), len(G2))
     1650#        print 'G: ', len(G),'G1: ', len(G1), 'G2: ', len(G2)
     1651
     1652        assert allclose(len(G), len(G1)+len(G2))
     1653        assert allclose(round(len(G)*factor), len(G1))
     1654       
     1655#        assert allclose(G == G1 + G2) must implentent __equal__
     1656       
    16291657         
    16301658if __name__ == "__main__":
    16311659
    1632     #suite = unittest.makeSuite(Test_Geospatial_data, 'test_ensure_geospatial')
     1660#    suite = unittest.makeSuite(Test_Geospatial_data, 'test_split')
     1661#    suite = unittest.makeSuite(Test_Geospatial_data, 'test_clip0')
    16331662    suite = unittest.makeSuite(Test_Geospatial_data, 'test')
    16341663    runner = unittest.TextTestRunner()
Note: See TracChangeset for help on using the changeset viewer.