source: branches/inundation-numpy-branch/geospatial_data/test_geospatial_data.py @ 6689

#!/usr/bin/env python


import unittest
import os
from numpy import zeros, array, allclose, concatenate
from math import sqrt, pi

from anuga.geospatial_data.geospatial_data import *

from coordinate_transforms.geo_reference import Geo_reference

class Test_Geospatial_data(unittest.TestCase):
    def setUp(self):
        pass

    def tearDown(self):
        pass


    def test_0(self):
        """Basic points
        """
        from coordinate_transforms.geo_reference import Geo_reference
        
        points = [[1.0, 2.1], [3.0, 5.3]]
        G = Geospatial_data(points)

        assert allclose(G.data_points, [[1.0, 2.1], [3.0, 5.3]])

        #Check defaults
        assert G.attributes is None
        
        assert G.geo_reference.zone == Geo_reference().zone
        assert G.geo_reference.xllcorner == Geo_reference().xllcorner
        assert G.geo_reference.yllcorner == Geo_reference().yllcorner
        

    def test_1(self):
        points = [[1.0, 2.1], [3.0, 5.3]]
        attributes = [2, 4]
        G = Geospatial_data(points, attributes)        

        assert G.attributes.keys()[0] == 'attribute'
        assert allclose(G.attributes.values()[0], [2, 4])
        

    def test_2(self):
        from coordinate_transforms.geo_reference import Geo_reference
        points = [[1.0, 2.1], [3.0, 5.3]]
        attributes = [2, 4]
        G = Geospatial_data(points, attributes,
                            geo_reference = Geo_reference(56, 100, 200))

        assert G.geo_reference.zone == 56
        assert G.geo_reference.xllcorner == 100
        assert G.geo_reference.yllcorner == 200


    def test_get_attributes_1(self):
        from coordinate_transforms.geo_reference import Geo_reference
        points = [[1.0, 2.1], [3.0, 5.3]]
        attributes = [2, 4]
        G = Geospatial_data(points, attributes,
                            geo_reference = Geo_reference(56, 100, 200))


        P = G.get_data_points(absolute = False)
        assert allclose(P, [[1.0, 2.1], [3.0, 5.3]])        

        P = G.get_data_points(absolute = True)
        assert allclose(P, [[101.0, 202.1], [103.0, 205.3]])        

        V = G.get_attributes() #Simply get them
        assert allclose(V, [2, 4])

        V = G.get_attributes('attribute') #Get by name
        assert allclose(V, [2, 4])

    def test_get_attributes_2(self):
        """Multiple attributes
        """
        
        from coordinate_transforms.geo_reference import Geo_reference
        points = [[1.0, 2.1], [3.0, 5.3]]
        attributes = {'a0': [0, 0], 'a1': [2, 4], 'a2': [79.4, -7]}
        G = Geospatial_data(points, attributes,
                            geo_reference = Geo_reference(56, 100, 200),
                            default_attribute_name = 'a1')


        P = G.get_data_points(absolute = False)
        assert allclose(P, [[1.0, 2.1], [3.0, 5.3]])        
        
        V = G.get_attributes() #Get default attribute
        assert allclose(V, [2, 4])

        V = G.get_attributes('a0') #Get by name
        assert allclose(V, [0, 0])

        V = G.get_attributes('a1') #Get by name
        assert allclose(V, [2, 4])

        V = G.get_attributes('a2') #Get by name
        assert allclose(V, [79.4, -7])

        try:
            V = G.get_attributes('hdnoatedu') #Invalid
        except AssertionError:
            pass
        else:
            raise 'Should have raised exception' 


    def test_conversions_to_points_dict(self):
        """test conversions to points_dict
        """
        
        from coordinate_transforms.geo_reference import Geo_reference
        points = [[1.0, 2.1], [3.0, 5.3]]
        attributes = {'a0': [0, 0], 'a1': [2, 4], 'a2': [79.4, -7]}
        G = Geospatial_data(points, attributes,
                            geo_reference = Geo_reference(56, 100, 200),
                            default_attribute_name = 'a1')


        points_dict = geospatial_data2points_dictionary(G)

        assert points_dict.has_key('pointlist')
        assert points_dict.has_key('attributelist')        
        assert points_dict.has_key('geo_reference')

        assert allclose( points_dict['pointlist'], points )

        A = points_dict['attributelist']
        assert A.has_key('a0')
        assert A.has_key('a1')
        assert A.has_key('a2')        

        assert allclose( A['a0'], [0, 0] )
        assert allclose( A['a1'], [2, 4] )        
        assert allclose( A['a2'], [79.4, -7] )


        geo = points_dict['geo_reference']
        assert geo is G.geo_reference


    def test_conversions_from_points_dict(self):
        """test conversions from points_dict
        """

        from coordinate_transforms.geo_reference import Geo_reference
        
        points = [[1.0, 2.1], [3.0, 5.3]]
        attributes = {'a0': [0, 0], 'a1': [2, 4], 'a2': [79.4, -7]}

        points_dict = {}
        points_dict['pointlist'] = points
        points_dict['attributelist'] = attributes
        points_dict['geo_reference'] = Geo_reference(56, 100, 200)
        

        G = points_dictionary2geospatial_data(points_dict)

        P = G.get_data_points(absolute = False)
        assert allclose(P, [[1.0, 2.1], [3.0, 5.3]])        
        
        #V = G.get_attribute_values() #Get default attribute
        #assert allclose(V, [2, 4])

        V = G.get_attributes('a0') #Get by name
        assert allclose(V, [0, 0])

        V = G.get_attributes('a1') #Get by name
        assert allclose(V, [2, 4])

        V = G.get_attributes('a2') #Get by name
        assert allclose(V, [79.4, -7])

    def test_add(self):
        """ test the addition of two geospatical objects
            no geo_reference see next test
        """
        points = [[1.0, 2.1], [3.0, 5.3]]
        attributes = {'depth':[2, 4], 'elevation':[6.1, 5]}
        attributes1 = {'depth':[2, 4], 'elevation':[2.5, 1]}
        G1 = Geospatial_data(points, attributes)        
        G2 = Geospatial_data(points, attributes1) 
        
#        print 'G1 depth=', G1.get_attributes('depth')
#        print 'G1 attrib keys=', G1.attributes.keys() 
        g3 = geospatial_data2points_dictionary(G1)
#        print 'g3=', g3
        
        G = G1 + G2

#        g = geospatial_data2points_dictionary(G)
#        print 'G=', g
#        print 'G points =', G.get_data_points()
#        print 'G attrib keys =', G.attributes.keys()
#        print 'G test =', G.get_attributes('elevation')
        assert G.attributes.has_key('depth')
        assert G.attributes.has_key('elevation')
        assert allclose(G.attributes['depth'], [2, 4, 2, 4])
        assert allclose(G.attributes['elevation'], [6.1, 5, 2.5, 1])
        assert allclose(G.get_data_points(), [[1.0, 2.1], [3.0, 5.3],
                                              [1.0, 2.1], [3.0, 5.3]])
        
    def test_add1 (self):
        """
        difference in Geo_reference resolved
        """
        points = [[1.0, 2.1], [3.0, 5.3]]
        points1 = [[5.0, 6.1], [6.0, 3.3]]
        attributes = [2, 4]
        attributes1 = [5, 76]
        geo_ref = Geo_reference(55, 1.0, 2.0)
        geo_ref1 = Geo_reference(zone = 55, xllcorner = 0.1,
                                yllcorner = 3.0, datum = 'wgs84',
                                projection = 'UTM', units = 'm')
                                
        G1 = Geospatial_data(points, attributes, geo_ref)
        G2 = Geospatial_data(points1, attributes1, geo_ref1)

        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]])                             

    def xtest_create_from_file(self):
        """Check that object can be created from a points file
        """

        from load_mesh.loadASCII import export_points_file
        

        points = [[1.0, 2.1], [3.0, 5.3], [5.0, 6.1], [6.0, 3.3]]
        attributes = [2, 4, 5, 76]

        # Use old pointsdict format
        pointsdict = {'pointlist': points,
                      'attributelist': {'att1': attributes,
                                        'att2': array(attributes) + 1}} 
        
        # Create points as an xya file 
        FN = 'test_points.xya'
        export_points_file(FN, pointsdict)


        #Create object from file
        G = Geospatial_data(filename = FN)

        assert allclose(G.get_data_points(), points)
        assert allclose(G.get_attributes('att1'), attributes)
        assert allclose(G.get_attributes('att2'), array(attributes) + 1)
        

        
        os.remove(FN)
        

if __name__ == "__main__":
    suite = unittest.makeSuite(Test_Geospatial_data,'test')
    runner = unittest.TextTestRunner()
    runner.run(suite)