#!/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)