[6086] | 1 | #!/usr/bin/env python |
---|
| 2 | |
---|
| 3 | import unittest |
---|
| 4 | import os |
---|
[6304] | 5 | import tempfile |
---|
[6086] | 6 | from math import sqrt, pi |
---|
| 7 | from sets import ImmutableSet |
---|
| 8 | |
---|
[6304] | 9 | import numpy as num |
---|
[6153] | 10 | |
---|
[6086] | 11 | from anuga.geospatial_data.geospatial_data import * |
---|
| 12 | from anuga.coordinate_transforms.geo_reference import Geo_reference, TitleError |
---|
| 13 | from anuga.coordinate_transforms.redfearn import degminsec2decimal_degrees |
---|
| 14 | from anuga.utilities.anuga_exceptions import ANUGAError |
---|
| 15 | from anuga.utilities.system_tools import get_host_name |
---|
| 16 | from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a |
---|
[6304] | 17 | from anuga.config import netcdf_float |
---|
[6086] | 18 | |
---|
[6304] | 19 | |
---|
[6086] | 20 | class Test_Geospatial_data(unittest.TestCase): |
---|
| 21 | def setUp(self): |
---|
| 22 | pass |
---|
| 23 | |
---|
| 24 | def tearDown(self): |
---|
| 25 | pass |
---|
| 26 | |
---|
| 27 | def test_0(self): |
---|
| 28 | #Basic points |
---|
| 29 | from anuga.coordinate_transforms.geo_reference import Geo_reference |
---|
[6304] | 30 | |
---|
[6086] | 31 | points = [[1.0, 2.1], [3.0, 5.3]] |
---|
| 32 | G = Geospatial_data(points) |
---|
[6153] | 33 | assert num.allclose(G.data_points, [[1.0, 2.1], [3.0, 5.3]]) |
---|
[6086] | 34 | |
---|
| 35 | # Check __repr__ |
---|
| 36 | # FIXME (Ole): Is this really machine independent? |
---|
| 37 | rep = `G` |
---|
| 38 | ref = '[[ 1. 2.1]\n [ 3. 5.3]]' |
---|
[6304] | 39 | msg = 'Representation %s is not equal to %s' % (rep, ref) |
---|
[6086] | 40 | assert rep == ref, msg |
---|
| 41 | |
---|
| 42 | #Check getter |
---|
[6153] | 43 | assert num.allclose(G.get_data_points(), [[1.0, 2.1], [3.0, 5.3]]) |
---|
[6304] | 44 | |
---|
[6086] | 45 | #Check defaults |
---|
| 46 | assert G.attributes is None |
---|
| 47 | assert G.geo_reference.zone == Geo_reference().zone |
---|
| 48 | assert G.geo_reference.xllcorner == Geo_reference().xllcorner |
---|
| 49 | assert G.geo_reference.yllcorner == Geo_reference().yllcorner |
---|
| 50 | |
---|
| 51 | def test_1(self): |
---|
| 52 | points = [[1.0, 2.1], [3.0, 5.3]] |
---|
| 53 | attributes = [2, 4] |
---|
[6304] | 54 | G = Geospatial_data(points, attributes) |
---|
[6086] | 55 | assert G.attributes.keys()[0] == DEFAULT_ATTRIBUTE |
---|
[6153] | 56 | assert num.allclose(G.attributes.values()[0], [2, 4]) |
---|
[6086] | 57 | |
---|
| 58 | def test_2(self): |
---|
| 59 | from anuga.coordinate_transforms.geo_reference import Geo_reference |
---|
[6304] | 60 | |
---|
[6086] | 61 | points = [[1.0, 2.1], [3.0, 5.3]] |
---|
| 62 | attributes = [2, 4] |
---|
| 63 | G = Geospatial_data(points, attributes, |
---|
| 64 | geo_reference=Geo_reference(56, 100, 200)) |
---|
| 65 | |
---|
| 66 | assert G.geo_reference.zone == 56 |
---|
| 67 | assert G.geo_reference.xllcorner == 100 |
---|
| 68 | assert G.geo_reference.yllcorner == 200 |
---|
| 69 | |
---|
| 70 | def test_get_attributes_1(self): |
---|
| 71 | from anuga.coordinate_transforms.geo_reference import Geo_reference |
---|
[6304] | 72 | |
---|
[6086] | 73 | points = [[1.0, 2.1], [3.0, 5.3]] |
---|
| 74 | attributes = [2, 4] |
---|
| 75 | G = Geospatial_data(points, attributes, |
---|
| 76 | geo_reference=Geo_reference(56, 100, 200)) |
---|
| 77 | |
---|
| 78 | P = G.get_data_points(absolute=False) |
---|
[6304] | 79 | assert num.allclose(P, [[1.0, 2.1], [3.0, 5.3]]) |
---|
[6086] | 80 | |
---|
| 81 | P = G.get_data_points(absolute=True) |
---|
[6304] | 82 | assert num.allclose(P, [[101.0, 202.1], [103.0, 205.3]]) |
---|
[6086] | 83 | |
---|
| 84 | V = G.get_attributes() #Simply get them |
---|
[6153] | 85 | assert num.allclose(V, [2, 4]) |
---|
[6086] | 86 | |
---|
| 87 | V = G.get_attributes(DEFAULT_ATTRIBUTE) #Get by name |
---|
[6153] | 88 | assert num.allclose(V, [2, 4]) |
---|
[6086] | 89 | |
---|
| 90 | def test_get_attributes_2(self): |
---|
| 91 | #Multiple attributes |
---|
| 92 | from anuga.coordinate_transforms.geo_reference import Geo_reference |
---|
[6304] | 93 | |
---|
[6086] | 94 | points = [[1.0, 2.1], [3.0, 5.3]] |
---|
| 95 | attributes = {'a0': [0, 0], 'a1': [2, 4], 'a2': [79.4, -7]} |
---|
| 96 | G = Geospatial_data(points, attributes, |
---|
| 97 | geo_reference=Geo_reference(56, 100, 200), |
---|
| 98 | default_attribute_name='a1') |
---|
| 99 | |
---|
[6304] | 100 | P = G.get_data_points(absolute=False) |
---|
| 101 | assert num.allclose(P, [[1.0, 2.1], [3.0, 5.3]]) |
---|
[6086] | 102 | |
---|
| 103 | V = G.get_attributes() #Get default attribute |
---|
[6153] | 104 | assert num.allclose(V, [2, 4]) |
---|
[6086] | 105 | |
---|
| 106 | V = G.get_attributes('a0') #Get by name |
---|
[6153] | 107 | assert num.allclose(V, [0, 0]) |
---|
[6086] | 108 | |
---|
| 109 | V = G.get_attributes('a1') #Get by name |
---|
[6153] | 110 | assert num.allclose(V, [2, 4]) |
---|
[6086] | 111 | |
---|
| 112 | V = G.get_attributes('a2') #Get by name |
---|
[6153] | 113 | assert num.allclose(V, [79.4, -7]) |
---|
[6086] | 114 | |
---|
| 115 | try: |
---|
| 116 | V = G.get_attributes('hdnoatedu') #Invalid |
---|
| 117 | except AssertionError: |
---|
| 118 | pass |
---|
| 119 | else: |
---|
[6304] | 120 | raise Exception, 'Should have raised exception' |
---|
[6086] | 121 | |
---|
| 122 | def test_get_data_points(self): |
---|
[6304] | 123 | points_ab = [[12.5, 34.7], [-4.5, -60.0]] |
---|
[6086] | 124 | x_p = -10 |
---|
| 125 | y_p = -40 |
---|
| 126 | geo_ref = Geo_reference(56, x_p, y_p) |
---|
| 127 | points_rel = geo_ref.change_points_geo_ref(points_ab) |
---|
[6304] | 128 | |
---|
[6086] | 129 | spatial = Geospatial_data(points_rel, geo_reference=geo_ref) |
---|
| 130 | results = spatial.get_data_points(absolute=False) |
---|
[6153] | 131 | assert num.allclose(results, points_rel) |
---|
[6304] | 132 | |
---|
[6086] | 133 | x_p = -1770 |
---|
| 134 | y_p = 4.01 |
---|
| 135 | geo_ref = Geo_reference(56, x_p, y_p) |
---|
| 136 | points_rel = geo_ref.change_points_geo_ref(points_ab) |
---|
[6304] | 137 | results = spatial.get_data_points(geo_reference=geo_ref) |
---|
| 138 | |
---|
[6153] | 139 | assert num.allclose(results, points_rel) |
---|
[6086] | 140 | |
---|
| 141 | def test_get_data_points_lat_long(self): |
---|
[6304] | 142 | # lat long [-30.], [130] |
---|
| 143 | # Zone: 52 |
---|
| 144 | # Easting: 596450.153 Northing: 6680793.777 |
---|
| 145 | # lat long [-32.], [131] |
---|
| 146 | # Zone: 52 |
---|
| 147 | # Easting: 688927.638 Northing: 6457816.509 |
---|
[6086] | 148 | |
---|
[6304] | 149 | points_Lat_long = [[-30., 130], [-32, 131]] |
---|
| 150 | |
---|
| 151 | spatial = Geospatial_data(latitudes=[-30, -32.], longitudes=[130, 131]) |
---|
[6086] | 152 | results = spatial.get_data_points(as_lat_long=True) |
---|
[6153] | 153 | assert num.allclose(results, points_Lat_long) |
---|
[6304] | 154 | |
---|
[6086] | 155 | def test_get_data_points_lat_longII(self): |
---|
| 156 | # x,y North,east long,lat |
---|
| 157 | boundary_polygon = [[ 250000, 7630000]] |
---|
| 158 | zone = 50 |
---|
[6304] | 159 | |
---|
[6086] | 160 | geo_reference = Geo_reference(zone=zone) |
---|
[6304] | 161 | geo = Geospatial_data(boundary_polygon ,geo_reference=geo_reference) |
---|
[6086] | 162 | seg_lat_long = geo.get_data_points(as_lat_long=True) |
---|
[6304] | 163 | lat_result = degminsec2decimal_degrees(-21, 24, 54) |
---|
| 164 | long_result = degminsec2decimal_degrees(114, 35, 17.89) |
---|
| 165 | assert num.allclose(seg_lat_long[0][0], lat_result) #lat |
---|
| 166 | assert num.allclose(seg_lat_long[0][1], long_result) #long |
---|
[6086] | 167 | |
---|
| 168 | def test_get_data_points_lat_longIII(self): |
---|
| 169 | # x,y North,east long,lat |
---|
[6304] | 170 | # for northern hemisphere |
---|
[6086] | 171 | boundary_polygon = [[419944.8, 918642.4]] |
---|
| 172 | zone = 47 |
---|
[6304] | 173 | |
---|
[6086] | 174 | geo_reference = Geo_reference(zone=zone) |
---|
[6304] | 175 | geo = Geospatial_data(boundary_polygon, geo_reference=geo_reference) |
---|
[6086] | 176 | seg_lat_long = geo.get_data_points(as_lat_long=True, |
---|
| 177 | isSouthHemisphere=False) |
---|
[6304] | 178 | lat_result = degminsec2decimal_degrees(8.31, 0, 0) |
---|
| 179 | long_result = degminsec2decimal_degrees(98.273, 0, 0) |
---|
[6153] | 180 | assert num.allclose(seg_lat_long[0][0], lat_result)#lat |
---|
| 181 | assert num.allclose(seg_lat_long[0][1], long_result)#long |
---|
[6086] | 182 | |
---|
[6304] | 183 | def test_set_geo_reference(self): |
---|
| 184 | '''test_set_georeference |
---|
[6086] | 185 | |
---|
[6304] | 186 | Test that georeference can be changed without changing the |
---|
[6086] | 187 | absolute values. |
---|
[6304] | 188 | ''' |
---|
| 189 | |
---|
| 190 | points_ab = [[12.5, 34.7], [-4.5, -60.0]] |
---|
[6086] | 191 | x_p = -10 |
---|
| 192 | y_p = -40 |
---|
| 193 | geo_ref = Geo_reference(56, x_p, y_p) |
---|
| 194 | points_rel = geo_ref.change_points_geo_ref(points_ab) |
---|
[6304] | 195 | |
---|
[6086] | 196 | # Create without geo_ref properly set |
---|
[6304] | 197 | G = Geospatial_data(points_rel) |
---|
[6153] | 198 | assert not num.allclose(points_ab, G.get_data_points(absolute=True)) |
---|
[6304] | 199 | |
---|
[6086] | 200 | # Create the way it should be |
---|
| 201 | G = Geospatial_data(points_rel, geo_reference=geo_ref) |
---|
[6153] | 202 | assert num.allclose(points_ab, G.get_data_points(absolute=True)) |
---|
[6304] | 203 | |
---|
[6086] | 204 | # Change georeference and check that absolute values are unchanged. |
---|
| 205 | x_p = 10 |
---|
| 206 | y_p = 400 |
---|
| 207 | new_geo_ref = Geo_reference(56, x_p, y_p) |
---|
| 208 | G.set_geo_reference(new_geo_ref) |
---|
[6410] | 209 | msg = ('points_ab=\n%s\nbut G.get_data_points(absolute=True)=\n%s' |
---|
[6304] | 210 | % (str(points_ab), str(G.get_data_points(absolute=True)))) |
---|
| 211 | assert num.allclose(points_ab, G.get_data_points(absolute=True)), msg |
---|
[6086] | 212 | |
---|
| 213 | def test_conversions_to_points_dict(self): |
---|
| 214 | #test conversions to points_dict |
---|
| 215 | from anuga.coordinate_transforms.geo_reference import Geo_reference |
---|
[6304] | 216 | |
---|
[6086] | 217 | points = [[1.0, 2.1], [3.0, 5.3]] |
---|
| 218 | attributes = {'a0': [0, 0], 'a1': [2, 4], 'a2': [79.4, -7]} |
---|
| 219 | G = Geospatial_data(points, attributes, |
---|
| 220 | geo_reference=Geo_reference(56, 100, 200), |
---|
| 221 | default_attribute_name='a1') |
---|
| 222 | |
---|
| 223 | points_dict = geospatial_data2points_dictionary(G) |
---|
| 224 | |
---|
| 225 | assert points_dict.has_key('pointlist') |
---|
[6304] | 226 | assert points_dict.has_key('attributelist') |
---|
[6086] | 227 | assert points_dict.has_key('geo_reference') |
---|
| 228 | |
---|
[6153] | 229 | assert num.allclose( points_dict['pointlist'], points ) |
---|
[6086] | 230 | |
---|
| 231 | A = points_dict['attributelist'] |
---|
| 232 | assert A.has_key('a0') |
---|
| 233 | assert A.has_key('a1') |
---|
[6304] | 234 | assert A.has_key('a2') |
---|
[6086] | 235 | |
---|
[6153] | 236 | assert num.allclose( A['a0'], [0, 0] ) |
---|
[6304] | 237 | assert num.allclose( A['a1'], [2, 4] ) |
---|
[6153] | 238 | assert num.allclose( A['a2'], [79.4, -7] ) |
---|
[6086] | 239 | |
---|
| 240 | geo = points_dict['geo_reference'] |
---|
| 241 | assert geo is G.geo_reference |
---|
| 242 | |
---|
| 243 | def test_conversions_from_points_dict(self): |
---|
[6304] | 244 | '''test conversions from points_dict''' |
---|
[6086] | 245 | |
---|
| 246 | from anuga.coordinate_transforms.geo_reference import Geo_reference |
---|
[6304] | 247 | |
---|
[6086] | 248 | points = [[1.0, 2.1], [3.0, 5.3]] |
---|
| 249 | attributes = {'a0': [0, 0], 'a1': [2, 4], 'a2': [79.4, -7]} |
---|
| 250 | |
---|
| 251 | points_dict = {} |
---|
| 252 | points_dict['pointlist'] = points |
---|
| 253 | points_dict['attributelist'] = attributes |
---|
| 254 | points_dict['geo_reference'] = Geo_reference(56, 100, 200) |
---|
| 255 | |
---|
| 256 | G = points_dictionary2geospatial_data(points_dict) |
---|
| 257 | P = G.get_data_points(absolute=False) |
---|
[6304] | 258 | assert num.allclose(P, [[1.0, 2.1], [3.0, 5.3]]) |
---|
[6086] | 259 | |
---|
| 260 | V = G.get_attributes('a0') #Get by name |
---|
[6153] | 261 | assert num.allclose(V, [0, 0]) |
---|
[6086] | 262 | |
---|
| 263 | V = G.get_attributes('a1') #Get by name |
---|
[6153] | 264 | assert num.allclose(V, [2, 4]) |
---|
[6086] | 265 | |
---|
| 266 | V = G.get_attributes('a2') #Get by name |
---|
[6153] | 267 | assert num.allclose(V, [79.4, -7]) |
---|
[6086] | 268 | |
---|
| 269 | def test_add(self): |
---|
[6304] | 270 | '''test the addition of two geospatical objects |
---|
| 271 | no geo_reference see next test |
---|
| 272 | ''' |
---|
| 273 | |
---|
[6086] | 274 | points = [[1.0, 2.1], [3.0, 5.3]] |
---|
[6304] | 275 | attributes = {'depth': [2, 4], 'elevation': [6.1, 5]} |
---|
| 276 | attributes1 = {'depth': [2, 4], 'elevation': [2.5, 1]} |
---|
| 277 | G1 = Geospatial_data(points, attributes) |
---|
| 278 | G2 = Geospatial_data(points, attributes1) |
---|
| 279 | |
---|
[6086] | 280 | G = G1 + G2 |
---|
| 281 | |
---|
| 282 | assert G.attributes.has_key('depth') |
---|
| 283 | assert G.attributes.has_key('elevation') |
---|
[6153] | 284 | assert num.allclose(G.attributes['depth'], [2, 4, 2, 4]) |
---|
| 285 | assert num.allclose(G.attributes['elevation'], [6.1, 5, 2.5, 1]) |
---|
| 286 | assert num.allclose(G.get_data_points(), [[1.0, 2.1], [3.0, 5.3], |
---|
| 287 | [1.0, 2.1], [3.0, 5.3]]) |
---|
| 288 | |
---|
[6086] | 289 | def test_addII(self): |
---|
[6304] | 290 | '''test the addition of two geospatical objects |
---|
| 291 | no geo_reference see next test |
---|
| 292 | ''' |
---|
| 293 | |
---|
[6086] | 294 | points = [[1.0, 2.1], [3.0, 5.3]] |
---|
[6304] | 295 | attributes = {'depth': [2, 4]} |
---|
| 296 | G1 = Geospatial_data(points, attributes) |
---|
| 297 | |
---|
[6086] | 298 | points = [[5.0, 2.1], [3.0, 50.3]] |
---|
[6304] | 299 | attributes = {'depth': [200, 400]} |
---|
[6086] | 300 | G2 = Geospatial_data(points, attributes) |
---|
[6304] | 301 | |
---|
[6086] | 302 | G = G1 + G2 |
---|
| 303 | |
---|
[6304] | 304 | assert G.attributes.has_key('depth') |
---|
[6086] | 305 | assert G.attributes.keys(), ['depth'] |
---|
[6153] | 306 | assert num.allclose(G.attributes['depth'], [2, 4, 200, 400]) |
---|
| 307 | assert num.allclose(G.get_data_points(), [[1.0, 2.1], [3.0, 5.3], |
---|
| 308 | [5.0, 2.1], [3.0, 50.3]]) |
---|
[6304] | 309 | |
---|
[6086] | 310 | def test_add_with_geo (self): |
---|
[6304] | 311 | '''Difference in Geo_reference resolved''' |
---|
| 312 | |
---|
[6086] | 313 | points1 = [[1.0, 2.1], [3.0, 5.3]] |
---|
| 314 | points2 = [[5.0, 6.1], [6.0, 3.3]] |
---|
| 315 | attributes1 = [2, 4] |
---|
| 316 | attributes2 = [5, 76] |
---|
| 317 | geo_ref1= Geo_reference(55, 1.0, 2.0) |
---|
| 318 | geo_ref2 = Geo_reference(zone=55, |
---|
| 319 | xllcorner=0.1, |
---|
| 320 | yllcorner=3.0, |
---|
| 321 | datum='wgs84', |
---|
| 322 | projection='UTM', |
---|
| 323 | units='m') |
---|
[6304] | 324 | |
---|
[6086] | 325 | G1 = Geospatial_data(points1, attributes1, geo_ref1) |
---|
| 326 | G2 = Geospatial_data(points2, attributes2, geo_ref2) |
---|
| 327 | |
---|
| 328 | #Check that absolute values are as expected |
---|
| 329 | P1 = G1.get_data_points(absolute=True) |
---|
[6153] | 330 | assert num.allclose(P1, [[2.0, 4.1], [4.0, 7.3]]) |
---|
[6086] | 331 | |
---|
| 332 | P2 = G2.get_data_points(absolute=True) |
---|
[6304] | 333 | assert num.allclose(P2, [[5.1, 9.1], [6.1, 6.3]]) |
---|
| 334 | |
---|
[6086] | 335 | G = G1 + G2 |
---|
| 336 | |
---|
| 337 | # Check absoluteness |
---|
[6153] | 338 | assert num.allclose(G.get_geo_reference().get_xllcorner(), 0.0) |
---|
| 339 | assert num.allclose(G.get_geo_reference().get_yllcorner(), 0.0) |
---|
[6086] | 340 | |
---|
| 341 | P = G.get_data_points(absolute=True) |
---|
| 342 | |
---|
[6428] | 343 | assert num.allclose(P, num.concatenate((P1,P2), axis=0)) #??default# |
---|
[6410] | 344 | msg = ('P=\n%s\nshould be close to\n' |
---|
| 345 | '[[2.0, 4.1], [4.0, 7.3],\n' |
---|
| 346 | ' [5.1, 9.1], [6.1, 6.3]]' |
---|
| 347 | % str(P)) |
---|
[6153] | 348 | assert num.allclose(P, [[2.0, 4.1], [4.0, 7.3], |
---|
[6304] | 349 | [5.1, 9.1], [6.1, 6.3]]), msg |
---|
[6086] | 350 | |
---|
[6304] | 351 | def test_add_with_geo_absolute (self): |
---|
| 352 | '''Difference in Geo_reference resolved''' |
---|
[6086] | 353 | |
---|
[6153] | 354 | points1 = num.array([[2.0, 4.1], [4.0, 7.3]]) |
---|
[6304] | 355 | points2 = num.array([[5.1, 9.1], [6.1, 6.3]]) |
---|
[6086] | 356 | attributes1 = [2, 4] |
---|
| 357 | attributes2 = [5, 76] |
---|
| 358 | geo_ref1= Geo_reference(55, 1.0, 2.0) |
---|
| 359 | geo_ref2 = Geo_reference(55, 2.0, 3.0) |
---|
| 360 | |
---|
[6304] | 361 | G1 = Geospatial_data(points1 - [geo_ref1.get_xllcorner(), |
---|
| 362 | geo_ref1.get_yllcorner()], |
---|
[6086] | 363 | attributes1, geo_ref1) |
---|
[6304] | 364 | |
---|
| 365 | G2 = Geospatial_data(points2 - [geo_ref2.get_xllcorner(), |
---|
| 366 | geo_ref2.get_yllcorner()], |
---|
[6086] | 367 | attributes2, geo_ref2) |
---|
| 368 | |
---|
| 369 | #Check that absolute values are as expected |
---|
| 370 | P1 = G1.get_data_points(absolute=True) |
---|
[6153] | 371 | assert num.allclose(P1, points1) |
---|
[6086] | 372 | |
---|
| 373 | P1 = G1.get_data_points(absolute=False) |
---|
[6410] | 374 | msg = ('P1=\n%s\nbut points1 - <...>=\n%s' |
---|
[6304] | 375 | % (str(P1), |
---|
| 376 | str(points1 - [geo_ref1.get_xllcorner(), |
---|
| 377 | geo_ref1.get_yllcorner()]))) |
---|
| 378 | assert num.allclose(P1, points1 - [geo_ref1.get_xllcorner(), |
---|
| 379 | geo_ref1.get_yllcorner()]), msg |
---|
[6086] | 380 | |
---|
| 381 | P2 = G2.get_data_points(absolute=True) |
---|
[6153] | 382 | assert num.allclose(P2, points2) |
---|
[6086] | 383 | |
---|
| 384 | P2 = G2.get_data_points(absolute=False) |
---|
[6304] | 385 | assert num.allclose(P2, points2 - [geo_ref2.get_xllcorner(), |
---|
| 386 | geo_ref2.get_yllcorner()]) |
---|
| 387 | |
---|
[6086] | 388 | G = G1 + G2 |
---|
| 389 | P = G.get_data_points(absolute=True) |
---|
| 390 | |
---|
[6428] | 391 | assert num.allclose(P, num.concatenate((points1, points2), axis=0)) #??default# |
---|
[6086] | 392 | |
---|
[6304] | 393 | def test_add_with_None(self): |
---|
| 394 | '''test that None can be added to a geospatical objects''' |
---|
[6086] | 395 | |
---|
[6153] | 396 | points1 = num.array([[2.0, 4.1], [4.0, 7.3]]) |
---|
[6304] | 397 | points2 = num.array([[5.1, 9.1], [6.1, 6.3]]) |
---|
[6086] | 398 | |
---|
| 399 | geo_ref1= Geo_reference(55, 1.0, 2.0) |
---|
| 400 | geo_ref2 = Geo_reference(zone=55, |
---|
| 401 | xllcorner=0.1, |
---|
| 402 | yllcorner=3.0, |
---|
| 403 | datum='wgs84', |
---|
| 404 | projection='UTM', |
---|
| 405 | units='m') |
---|
| 406 | |
---|
[6304] | 407 | attributes1 = {'depth': [2, 4.7], 'elevation': [6.1, 5]} |
---|
| 408 | attributes2 = {'depth': [-2.3, 4], 'elevation': [2.5, 1]} |
---|
[6086] | 409 | |
---|
| 410 | G1 = Geospatial_data(points1, attributes1, geo_ref1) |
---|
[6153] | 411 | assert num.allclose(G1.get_geo_reference().get_xllcorner(), 1.0) |
---|
| 412 | assert num.allclose(G1.get_geo_reference().get_yllcorner(), 2.0) |
---|
[6086] | 413 | assert G1.attributes.has_key('depth') |
---|
| 414 | assert G1.attributes.has_key('elevation') |
---|
[6153] | 415 | assert num.allclose(G1.attributes['depth'], [2, 4.7]) |
---|
[6304] | 416 | assert num.allclose(G1.attributes['elevation'], [6.1, 5]) |
---|
| 417 | |
---|
[6086] | 418 | G2 = Geospatial_data(points2, attributes2, geo_ref2) |
---|
[6153] | 419 | assert num.allclose(G2.get_geo_reference().get_xllcorner(), 0.1) |
---|
| 420 | assert num.allclose(G2.get_geo_reference().get_yllcorner(), 3.0) |
---|
[6086] | 421 | assert G2.attributes.has_key('depth') |
---|
| 422 | assert G2.attributes.has_key('elevation') |
---|
[6153] | 423 | assert num.allclose(G2.attributes['depth'], [-2.3, 4]) |
---|
[6304] | 424 | assert num.allclose(G2.attributes['elevation'], [2.5, 1]) |
---|
[6086] | 425 | |
---|
| 426 | #Check that absolute values are as expected |
---|
| 427 | P1 = G1.get_data_points(absolute=True) |
---|
[6153] | 428 | assert num.allclose(P1, [[3.0, 6.1], [5.0, 9.3]]) |
---|
[6086] | 429 | |
---|
| 430 | P2 = G2.get_data_points(absolute=True) |
---|
[6304] | 431 | assert num.allclose(P2, [[5.2, 12.1], [6.2, 9.3]]) |
---|
[6086] | 432 | |
---|
| 433 | # Normal add |
---|
| 434 | G = G1 + None |
---|
| 435 | assert G.attributes.has_key('depth') |
---|
| 436 | assert G.attributes.has_key('elevation') |
---|
[6153] | 437 | assert num.allclose(G.attributes['depth'], [2, 4.7]) |
---|
[6304] | 438 | assert num.allclose(G.attributes['elevation'], [6.1, 5]) |
---|
[6086] | 439 | |
---|
| 440 | # Points are now absolute. |
---|
[6153] | 441 | assert num.allclose(G.get_geo_reference().get_xllcorner(), 0.0) |
---|
| 442 | assert num.allclose(G.get_geo_reference().get_yllcorner(), 0.0) |
---|
[6304] | 443 | P = G.get_data_points(absolute=True) |
---|
[6410] | 444 | msg = 'P=\n%s' % str(P) |
---|
[6304] | 445 | assert num.allclose(P, [[3.0, 6.1], [5.0, 9.3]]), msg |
---|
[6086] | 446 | |
---|
| 447 | G = G2 + None |
---|
| 448 | assert G.attributes.has_key('depth') |
---|
| 449 | assert G.attributes.has_key('elevation') |
---|
[6153] | 450 | assert num.allclose(G.attributes['depth'], [-2.3, 4]) |
---|
[6304] | 451 | assert num.allclose(G.attributes['elevation'], [2.5, 1]) |
---|
[6153] | 452 | assert num.allclose(G.get_geo_reference().get_xllcorner(), 0.0) |
---|
| 453 | assert num.allclose(G.get_geo_reference().get_yllcorner(), 0.0) |
---|
[6304] | 454 | |
---|
| 455 | P = G.get_data_points(absolute=True) |
---|
[6153] | 456 | assert num.allclose(P, [[5.2, 12.1], [6.2, 9.3]]) |
---|
[6086] | 457 | |
---|
| 458 | # Reverse add |
---|
| 459 | G = None + G1 |
---|
| 460 | assert G.attributes.has_key('depth') |
---|
| 461 | assert G.attributes.has_key('elevation') |
---|
[6153] | 462 | assert num.allclose(G.attributes['depth'], [2, 4.7]) |
---|
[6304] | 463 | assert num.allclose(G.attributes['elevation'], [6.1, 5]) |
---|
[6086] | 464 | |
---|
| 465 | # Points are now absolute. |
---|
[6153] | 466 | assert num.allclose(G.get_geo_reference().get_xllcorner(), 0.0) |
---|
| 467 | assert num.allclose(G.get_geo_reference().get_yllcorner(), 0.0) |
---|
[6086] | 468 | |
---|
[6304] | 469 | P = G.get_data_points(absolute=True) |
---|
| 470 | assert num.allclose(P, [[3.0, 6.1], [5.0, 9.3]]) |
---|
[6086] | 471 | |
---|
| 472 | G = None + G2 |
---|
| 473 | assert G.attributes.has_key('depth') |
---|
| 474 | assert G.attributes.has_key('elevation') |
---|
[6153] | 475 | assert num.allclose(G.attributes['depth'], [-2.3, 4]) |
---|
[6304] | 476 | assert num.allclose(G.attributes['elevation'], [2.5, 1]) |
---|
[6086] | 477 | |
---|
[6153] | 478 | assert num.allclose(G.get_geo_reference().get_xllcorner(), 0.0) |
---|
| 479 | assert num.allclose(G.get_geo_reference().get_yllcorner(), 0.0) |
---|
[6304] | 480 | |
---|
| 481 | P = G.get_data_points(absolute=True) |
---|
[6153] | 482 | assert num.allclose(P, [[5.2, 12.1], [6.2, 9.3]]) |
---|
[6086] | 483 | |
---|
[6304] | 484 | def test_clip0(self): |
---|
| 485 | '''test_clip0(self): |
---|
[6086] | 486 | |
---|
| 487 | Test that point sets can be clipped by a polygon |
---|
[6304] | 488 | ''' |
---|
| 489 | |
---|
[6086] | 490 | from anuga.coordinate_transforms.geo_reference import Geo_reference |
---|
[6304] | 491 | |
---|
[6086] | 492 | points = [[-1, 4], [0.2, 0.5], [1.0, 2.1], [0.4, 0.3], [3.0, 5.3], |
---|
| 493 | [0, 0], [2.4, 3.3]] |
---|
| 494 | G = Geospatial_data(points) |
---|
| 495 | |
---|
[6304] | 496 | # First try the unit square |
---|
| 497 | U = [[0,0], [1,0], [1,1], [0,1]] |
---|
| 498 | assert num.allclose(G.clip(U).get_data_points(), |
---|
| 499 | [[0.2, 0.5], [0.4, 0.3], [0, 0]]) |
---|
[6086] | 500 | |
---|
| 501 | # Then a more complex polygon |
---|
| 502 | polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]] |
---|
[6304] | 503 | points = [[0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], |
---|
| 504 | [0.5, 1.5], [0.5, -0.5]] |
---|
[6086] | 505 | G = Geospatial_data(points) |
---|
| 506 | |
---|
[6153] | 507 | assert num.allclose(G.clip(polygon).get_data_points(), |
---|
[6304] | 508 | [[0.5, 0.5], [1, -0.5], [1.5, 0]]) |
---|
[6086] | 509 | |
---|
| 510 | def test_clip0_with_attributes(self): |
---|
[6304] | 511 | '''test_clip0_with_attributes(self): |
---|
| 512 | |
---|
[6086] | 513 | Test that point sets with attributes can be clipped by a polygon |
---|
[6304] | 514 | ''' |
---|
| 515 | |
---|
[6086] | 516 | from anuga.coordinate_transforms.geo_reference import Geo_reference |
---|
[6304] | 517 | |
---|
[6086] | 518 | points = [[-1, 4], [0.2, 0.5], [1.0, 2.1], [0.4, 0.3], [3.0, 5.3], |
---|
| 519 | [0, 0], [2.4, 3.3]] |
---|
| 520 | |
---|
| 521 | attributes = [2, -4, 5, 76, -2, 0.1, 3] |
---|
| 522 | att_dict = {'att1': attributes, |
---|
[6153] | 523 | 'att2': num.array(attributes)+1} |
---|
[6304] | 524 | |
---|
[6086] | 525 | G = Geospatial_data(points, att_dict) |
---|
| 526 | |
---|
[6304] | 527 | # First try the unit square |
---|
| 528 | U = [[0,0], [1,0], [1,1], [0,1]] |
---|
| 529 | assert num.allclose(G.clip(U).get_data_points(), |
---|
| 530 | [[0.2, 0.5], [0.4, 0.3], [0, 0]]) |
---|
[6153] | 531 | assert num.allclose(G.clip(U).get_attributes('att1'), [-4, 76, 0.1]) |
---|
[6304] | 532 | assert num.allclose(G.clip(U).get_attributes('att2'), [-3, 77, 1.1]) |
---|
[6086] | 533 | |
---|
| 534 | # Then a more complex polygon |
---|
| 535 | polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]] |
---|
[6304] | 536 | points = [[0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], |
---|
| 537 | [0.5, 1.5], [0.5, -0.5]] |
---|
[6086] | 538 | |
---|
| 539 | # This time just one attribute |
---|
| 540 | attributes = [2, -4, 5, 76, -2, 0.1] |
---|
| 541 | G = Geospatial_data(points, attributes) |
---|
[6153] | 542 | assert num.allclose(G.clip(polygon).get_data_points(), |
---|
[6304] | 543 | [[0.5, 0.5], [1, -0.5], [1.5, 0]]) |
---|
[6153] | 544 | assert num.allclose(G.clip(polygon).get_attributes(), [-4, 5, 76]) |
---|
[6086] | 545 | |
---|
| 546 | def test_clip1(self): |
---|
[6304] | 547 | '''test_clip1(self): |
---|
| 548 | |
---|
[6086] | 549 | Test that point sets can be clipped by a polygon given as |
---|
| 550 | another Geospatial dataset |
---|
[6304] | 551 | ''' |
---|
| 552 | |
---|
[6086] | 553 | from anuga.coordinate_transforms.geo_reference import Geo_reference |
---|
[6304] | 554 | |
---|
[6086] | 555 | points = [[-1, 4], [0.2, 0.5], [1.0, 2.1], [0.4, 0.3], [3.0, 5.3], |
---|
| 556 | [0, 0], [2.4, 3.3]] |
---|
| 557 | attributes = [2, -4, 5, 76, -2, 0.1, 3] |
---|
| 558 | att_dict = {'att1': attributes, |
---|
[6153] | 559 | 'att2': num.array(attributes)+1} |
---|
[6086] | 560 | G = Geospatial_data(points, att_dict) |
---|
[6304] | 561 | |
---|
| 562 | # First try the unit square |
---|
| 563 | U = Geospatial_data([[0,0], [1,0], [1,1], [0,1]]) |
---|
[6153] | 564 | assert num.allclose(G.clip(U).get_data_points(), |
---|
[6304] | 565 | [[0.2, 0.5], [0.4, 0.3], [0, 0]]) |
---|
| 566 | assert num.allclose(G.clip(U).get_attributes('att1'), [-4, 76, 0.1]) |
---|
| 567 | assert num.allclose(G.clip(U).get_attributes('att2'), [-3, 77, 1.1]) |
---|
[6086] | 568 | |
---|
| 569 | # Then a more complex polygon |
---|
[6304] | 570 | points = [[0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], |
---|
| 571 | [0.5, 1.5], [0.5, -0.5]] |
---|
| 572 | attributes = [2, -4, 5, 76, -2, 0.1] |
---|
[6086] | 573 | G = Geospatial_data(points, attributes) |
---|
[6304] | 574 | polygon = Geospatial_data([[0,0], [1,0], [0.5,-1], [2, -1], |
---|
| 575 | [2,1], [0,1]]) |
---|
[6086] | 576 | |
---|
[6153] | 577 | assert num.allclose(G.clip(polygon).get_data_points(), |
---|
| 578 | [[0.5, 0.5], [1, -0.5], [1.5, 0]]) |
---|
| 579 | assert num.allclose(G.clip(polygon).get_attributes(), [-4, 5, 76]) |
---|
[6086] | 580 | |
---|
| 581 | def test_clip0_outside(self): |
---|
[6304] | 582 | '''test_clip0_outside(self): |
---|
| 583 | |
---|
[6086] | 584 | Test that point sets can be clipped outside of a polygon |
---|
[6304] | 585 | ''' |
---|
| 586 | |
---|
[6086] | 587 | from anuga.coordinate_transforms.geo_reference import Geo_reference |
---|
[6304] | 588 | |
---|
[6086] | 589 | points = [[-1, 4], [0.2, 0.5], [1.0, 2.1], [0.4, 0.3], [3.0, 5.3], |
---|
| 590 | [0, 0], [2.4, 3.3]] |
---|
[6304] | 591 | attributes = [2, -4, 5, 76, -2, 0.1, 3] |
---|
[6086] | 592 | G = Geospatial_data(points, attributes) |
---|
| 593 | |
---|
[6304] | 594 | # First try the unit square |
---|
[6086] | 595 | U = [[0,0], [1,0], [1,1], [0,1]] |
---|
[6153] | 596 | assert num.allclose(G.clip_outside(U).get_data_points(), |
---|
| 597 | [[-1, 4], [1.0, 2.1], [3.0, 5.3], [2.4, 3.3]]) |
---|
[6304] | 598 | assert num.allclose(G.clip_outside(U).get_attributes(), [2, 5, -2, 3]) |
---|
[6086] | 599 | |
---|
| 600 | # Then a more complex polygon |
---|
| 601 | polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]] |
---|
[6304] | 602 | points = [[0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], |
---|
| 603 | [0.5, 1.5], [0.5, -0.5]] |
---|
| 604 | attributes = [2, -4, 5, 76, -2, 0.1] |
---|
[6086] | 605 | G = Geospatial_data(points, attributes) |
---|
[6153] | 606 | assert num.allclose(G.clip_outside(polygon).get_data_points(), |
---|
| 607 | [[0.5, 1.4], [0.5, 1.5], [0.5, -0.5]]) |
---|
[6304] | 608 | assert num.allclose(G.clip_outside(polygon).get_attributes(), |
---|
| 609 | [2, -2, 0.1]) |
---|
[6086] | 610 | |
---|
[6304] | 611 | def test_clip1_outside(self): |
---|
| 612 | '''test_clip1_outside(self): |
---|
[6086] | 613 | |
---|
| 614 | Test that point sets can be clipped outside of a polygon given as |
---|
| 615 | another Geospatial dataset |
---|
[6304] | 616 | ''' |
---|
| 617 | |
---|
[6086] | 618 | from anuga.coordinate_transforms.geo_reference import Geo_reference |
---|
[6304] | 619 | |
---|
[6086] | 620 | points = [[-1, 4], [0.2, 0.5], [1.0, 2.1], [0.4, 0.3], [3.0, 5.3], |
---|
| 621 | [0, 0], [2.4, 3.3]] |
---|
[6304] | 622 | attributes = [2, -4, 5, 76, -2, 0.1, 3] |
---|
| 623 | G = Geospatial_data(points, attributes) |
---|
[6086] | 624 | |
---|
[6304] | 625 | # First try the unit square |
---|
| 626 | U = Geospatial_data([[0,0], [1,0], [1,1], [0,1]]) |
---|
[6153] | 627 | assert num.allclose(G.clip_outside(U).get_data_points(), |
---|
| 628 | [[-1, 4], [1.0, 2.1], [3.0, 5.3], [2.4, 3.3]]) |
---|
[6304] | 629 | assert num.allclose(G.clip(U).get_attributes(), [-4, 76, 0.1]) |
---|
[6086] | 630 | |
---|
| 631 | # Then a more complex polygon |
---|
[6304] | 632 | points = [[0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], |
---|
| 633 | [0.5, 1.5], [0.5, -0.5]] |
---|
| 634 | attributes = [2, -4, 5, 76, -2, 0.1] |
---|
[6086] | 635 | G = Geospatial_data(points, attributes) |
---|
[6304] | 636 | polygon = Geospatial_data([[0, 0], [1, 0], [0.5, -1], [2, -1], |
---|
| 637 | [2, 1], [0, 1]]) |
---|
[6153] | 638 | assert num.allclose(G.clip_outside(polygon).get_data_points(), |
---|
| 639 | [[0.5, 1.4], [0.5, 1.5], [0.5, -0.5]]) |
---|
[6304] | 640 | assert num.allclose(G.clip_outside(polygon).get_attributes(), |
---|
| 641 | [2, -2, 0.1]) |
---|
[6086] | 642 | |
---|
[6304] | 643 | def test_clip1_inside_outside(self): |
---|
| 644 | '''test_clip1_inside_outside(self): |
---|
[6086] | 645 | |
---|
| 646 | Test that point sets can be clipped outside of a polygon given as |
---|
| 647 | another Geospatial dataset |
---|
[6304] | 648 | ''' |
---|
| 649 | |
---|
[6086] | 650 | from anuga.coordinate_transforms.geo_reference import Geo_reference |
---|
[6304] | 651 | |
---|
[6086] | 652 | points = [[-1, 4], [0.2, 0.5], [1.0, 2.1], [0.4, 0.3], [3.0, 5.3], |
---|
| 653 | [0, 0], [2.4, 3.3]] |
---|
[6304] | 654 | attributes = [2, -4, 5, 76, -2, 0.1, 3] |
---|
[6086] | 655 | G = Geospatial_data(points, attributes) |
---|
| 656 | |
---|
[6304] | 657 | # First try the unit square |
---|
| 658 | U = Geospatial_data([[0,0], [1,0], [1,1], [0,1]]) |
---|
[6086] | 659 | G1 = G.clip(U) |
---|
[6304] | 660 | assert num.allclose(G1.get_data_points(), |
---|
| 661 | [[0.2, 0.5], [0.4, 0.3], [0, 0]]) |
---|
[6153] | 662 | assert num.allclose(G.clip(U).get_attributes(), [-4, 76, 0.1]) |
---|
[6086] | 663 | G2 = G.clip_outside(U) |
---|
[6153] | 664 | assert num.allclose(G2.get_data_points(),[[-1, 4], [1.0, 2.1], |
---|
| 665 | [3.0, 5.3], [2.4, 3.3]]) |
---|
[6304] | 666 | assert num.allclose(G.clip_outside(U).get_attributes(), [2, 5, -2, 3]) |
---|
[6086] | 667 | |
---|
| 668 | # New ordering |
---|
[6304] | 669 | new_points = [[0.2, 0.5], [0.4, 0.3], [0, 0], [-1, 4], |
---|
| 670 | [1.0, 2.1], [3.0, 5.3], [2.4, 3.3]] |
---|
| 671 | new_attributes = [-4, 76, 0.1, 2, 5, -2, 3] |
---|
[6086] | 672 | |
---|
[6153] | 673 | assert num.allclose((G1+G2).get_data_points(), new_points) |
---|
| 674 | assert num.allclose((G1+G2).get_attributes(), new_attributes) |
---|
[6086] | 675 | |
---|
| 676 | G = G1+G2 |
---|
| 677 | FN = 'test_combine.pts' |
---|
| 678 | G.export_points_file(FN) |
---|
| 679 | |
---|
| 680 | # Read it back in |
---|
| 681 | G3 = Geospatial_data(FN) |
---|
| 682 | |
---|
[6304] | 683 | # Check result |
---|
| 684 | assert num.allclose(G3.get_data_points(), new_points) |
---|
| 685 | assert num.allclose(G3.get_attributes(), new_attributes) |
---|
[6086] | 686 | |
---|
| 687 | os.remove(FN) |
---|
| 688 | |
---|
| 689 | def test_load_csv(self): |
---|
[6304] | 690 | fileName = tempfile.mktemp('.csv') |
---|
| 691 | file = open(fileName,'w') |
---|
| 692 | file.write('x,y,elevation speed \n\ |
---|
[6086] | 693 | 1.0 0.0 10.0 0.0\n\ |
---|
| 694 | 0.0 1.0 0.0 10.0\n\ |
---|
[6304] | 695 | 1.0 0.0 10.4 40.0\n') |
---|
[6086] | 696 | file.close() |
---|
[6304] | 697 | |
---|
[6086] | 698 | results = Geospatial_data(fileName) |
---|
| 699 | os.remove(fileName) |
---|
[6304] | 700 | assert num.allclose(results.get_data_points(), |
---|
| 701 | [[1.0, 0.0],[0.0, 1.0], [1.0, 0.0]]) |
---|
[6153] | 702 | assert num.allclose(results.get_attributes(attribute_name='elevation'), |
---|
| 703 | [10.0, 0.0, 10.4]) |
---|
| 704 | assert num.allclose(results.get_attributes(attribute_name='speed'), |
---|
| 705 | [0.0, 10.0, 40.0]) |
---|
[6086] | 706 | |
---|
[6304] | 707 | ################################################################################ |
---|
| 708 | # Test CSV files |
---|
| 709 | ################################################################################ |
---|
[6086] | 710 | |
---|
| 711 | def test_load_csv_lat_long_bad_blocking(self): |
---|
[6304] | 712 | '''test_load_csv_lat_long_bad_blocking(self): |
---|
[6086] | 713 | Different zones in Geo references |
---|
[6304] | 714 | ''' |
---|
| 715 | |
---|
| 716 | fileName = tempfile.mktemp('.csv') |
---|
| 717 | file = open(fileName, 'w') |
---|
| 718 | file.write('Lati,LONG,z \n\ |
---|
[6086] | 719 | -25.0,180.0,452.688000\n\ |
---|
[6304] | 720 | -34,150.0,459.126000\n') |
---|
[6086] | 721 | file.close() |
---|
[6304] | 722 | |
---|
[6086] | 723 | results = Geospatial_data(fileName, max_read_lines=1, |
---|
| 724 | load_file_now=False) |
---|
[6304] | 725 | |
---|
[6086] | 726 | try: |
---|
| 727 | for i in results: |
---|
| 728 | pass |
---|
| 729 | except ANUGAError: |
---|
| 730 | pass |
---|
| 731 | else: |
---|
| 732 | msg = 'Different zones in Geo references not caught.' |
---|
[6304] | 733 | raise Exception, msg |
---|
| 734 | |
---|
[6086] | 735 | os.remove(fileName) |
---|
[6304] | 736 | |
---|
[6086] | 737 | def test_load_csv(self): |
---|
[6304] | 738 | fileName = tempfile.mktemp('.txt') |
---|
| 739 | file = open(fileName, 'w') |
---|
| 740 | file.write(' x,y, elevation , speed \n\ |
---|
[6086] | 741 | 1.0, 0.0, 10.0, 0.0\n\ |
---|
| 742 | 0.0, 1.0, 0.0, 10.0\n\ |
---|
[6304] | 743 | 1.0, 0.0 ,10.4, 40.0\n') |
---|
[6086] | 744 | file.close() |
---|
| 745 | |
---|
| 746 | results = Geospatial_data(fileName, max_read_lines=2) |
---|
| 747 | |
---|
[6304] | 748 | assert num.allclose(results.get_data_points(), |
---|
| 749 | [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]]) |
---|
| 750 | assert num.allclose(results.get_attributes(attribute_name='elevation'), |
---|
| 751 | [10.0, 0.0, 10.4]) |
---|
| 752 | assert num.allclose(results.get_attributes(attribute_name='speed'), |
---|
| 753 | [0.0, 10.0, 40.0]) |
---|
[6086] | 754 | |
---|
| 755 | # Blocking |
---|
| 756 | geo_list = [] |
---|
| 757 | for i in results: |
---|
| 758 | geo_list.append(i) |
---|
[6304] | 759 | |
---|
[6153] | 760 | assert num.allclose(geo_list[0].get_data_points(), |
---|
| 761 | [[1.0, 0.0],[0.0, 1.0]]) |
---|
| 762 | assert num.allclose(geo_list[0].get_attributes(attribute_name='elevation'), |
---|
| 763 | [10.0, 0.0]) |
---|
| 764 | assert num.allclose(geo_list[1].get_data_points(), |
---|
[6304] | 765 | [[1.0, 0.0]]) |
---|
[6153] | 766 | assert num.allclose(geo_list[1].get_attributes(attribute_name='elevation'), |
---|
| 767 | [10.4]) |
---|
[6304] | 768 | |
---|
| 769 | os.remove(fileName) |
---|
| 770 | |
---|
[6086] | 771 | def test_load_csv_bad(self): |
---|
[6304] | 772 | '''test_load_csv_bad(self): |
---|
[6086] | 773 | header column, body column missmatch |
---|
| 774 | (Format error) |
---|
[6304] | 775 | ''' |
---|
| 776 | |
---|
| 777 | fileName = tempfile.mktemp('.txt') |
---|
| 778 | file = open(fileName, 'w') |
---|
| 779 | file.write(' elevation , speed \n\ |
---|
[6086] | 780 | 1.0, 0.0, 10.0, 0.0\n\ |
---|
| 781 | 0.0, 1.0, 0.0, 10.0\n\ |
---|
[6304] | 782 | 1.0, 0.0 ,10.4, 40.0\n') |
---|
[6086] | 783 | file.close() |
---|
| 784 | |
---|
| 785 | results = Geospatial_data(fileName, max_read_lines=2, |
---|
| 786 | load_file_now=False) |
---|
| 787 | |
---|
| 788 | # Blocking |
---|
| 789 | geo_list = [] |
---|
| 790 | try: |
---|
| 791 | for i in results: |
---|
| 792 | geo_list.append(i) |
---|
| 793 | except SyntaxError: |
---|
| 794 | pass |
---|
| 795 | else: |
---|
[6304] | 796 | msg = 'Bad file did not raise error!' |
---|
| 797 | raise Exception, msg |
---|
[6086] | 798 | os.remove(fileName) |
---|
| 799 | |
---|
| 800 | def test_load_csv_badII(self): |
---|
[6304] | 801 | '''test_load_csv_bad(self): |
---|
[6086] | 802 | header column, body column missmatch |
---|
| 803 | (Format error) |
---|
[6304] | 804 | ''' |
---|
| 805 | |
---|
| 806 | fileName = tempfile.mktemp('.txt') |
---|
| 807 | file = open(fileName, 'w') |
---|
| 808 | file.write(' x,y,elevation , speed \n\ |
---|
[6086] | 809 | 1.0, 0.0, 10.0, 0.0\n\ |
---|
| 810 | 0.0, 1.0, 0.0, 10\n\ |
---|
[6304] | 811 | 1.0, 0.0 ,10.4 yeah\n') |
---|
[6086] | 812 | file.close() |
---|
| 813 | |
---|
| 814 | results = Geospatial_data(fileName, max_read_lines=2, |
---|
| 815 | load_file_now=False) |
---|
| 816 | |
---|
| 817 | # Blocking |
---|
| 818 | geo_list = [] |
---|
| 819 | try: |
---|
| 820 | for i in results: |
---|
| 821 | geo_list.append(i) |
---|
| 822 | except SyntaxError: |
---|
| 823 | pass |
---|
| 824 | else: |
---|
[6304] | 825 | msg = 'Bad file did not raise error!' |
---|
| 826 | raise Exception, msg |
---|
[6086] | 827 | os.remove(fileName) |
---|
| 828 | |
---|
| 829 | def test_load_csv_badIII(self): |
---|
[6304] | 830 | '''test_load_csv_bad(self): |
---|
[6086] | 831 | space delimited |
---|
[6304] | 832 | ''' |
---|
| 833 | |
---|
| 834 | fileName = tempfile.mktemp('.txt') |
---|
| 835 | file = open(fileName, 'w') |
---|
| 836 | file.write(' x y elevation speed \n\ |
---|
[6086] | 837 | 1. 0.0 10.0 0.0\n\ |
---|
| 838 | 0.0 1.0 0.0 10.0\n\ |
---|
[6304] | 839 | 1.0 0.0 10.4 40.0\n') |
---|
[6086] | 840 | file.close() |
---|
| 841 | |
---|
| 842 | try: |
---|
| 843 | results = Geospatial_data(fileName, max_read_lines=2, |
---|
| 844 | load_file_now=True) |
---|
| 845 | except SyntaxError: |
---|
| 846 | pass |
---|
| 847 | else: |
---|
[6304] | 848 | msg = 'Bad file did not raise error!' |
---|
| 849 | raise Exception, msg |
---|
[6086] | 850 | os.remove(fileName) |
---|
[6304] | 851 | |
---|
[6086] | 852 | def test_load_csv_badIV(self): |
---|
[6304] | 853 | ''' test_load_csv_bad(self): |
---|
[6086] | 854 | header column, body column missmatch |
---|
| 855 | (Format error) |
---|
[6304] | 856 | ''' |
---|
| 857 | |
---|
| 858 | fileName = tempfile.mktemp('.txt') |
---|
| 859 | file = open(fileName, 'w') |
---|
| 860 | file.write(' x,y,elevation , speed \n\ |
---|
[6086] | 861 | 1.0, 0.0, 10.0, wow\n\ |
---|
| 862 | 0.0, 1.0, 0.0, ha\n\ |
---|
[6304] | 863 | 1.0, 0.0 ,10.4, yeah\n') |
---|
[6086] | 864 | file.close() |
---|
| 865 | |
---|
| 866 | results = Geospatial_data(fileName, max_read_lines=2, |
---|
| 867 | load_file_now=False) |
---|
| 868 | |
---|
| 869 | # Blocking |
---|
| 870 | geo_list = [] |
---|
| 871 | try: |
---|
| 872 | for i in results: |
---|
| 873 | geo_list.append(i) |
---|
| 874 | except SyntaxError: |
---|
| 875 | pass |
---|
| 876 | else: |
---|
[6304] | 877 | msg = 'Bad file did not raise error!' |
---|
| 878 | raise Exception, msg |
---|
[6086] | 879 | os.remove(fileName) |
---|
| 880 | |
---|
| 881 | def test_load_pts_blocking(self): |
---|
| 882 | #This is pts! |
---|
[6304] | 883 | fileName = tempfile.mktemp('.txt') |
---|
| 884 | file = open(fileName, 'w') |
---|
| 885 | file.write(' x,y, elevation , speed \n\ |
---|
[6086] | 886 | 1.0, 0.0, 10.0, 0.0\n\ |
---|
| 887 | 0.0, 1.0, 0.0, 10.0\n\ |
---|
[6304] | 888 | 1.0, 0.0 ,10.4, 40.0\n') |
---|
[6086] | 889 | file.close() |
---|
| 890 | |
---|
[6304] | 891 | pts_file = tempfile.mktemp('.pts') |
---|
| 892 | |
---|
[6086] | 893 | convert = Geospatial_data(fileName) |
---|
| 894 | convert.export_points_file(pts_file) |
---|
| 895 | results = Geospatial_data(pts_file, max_read_lines=2) |
---|
| 896 | |
---|
[6304] | 897 | assert num.allclose(results.get_data_points(), |
---|
| 898 | [[1.0, 0.0],[0.0, 1.0], [1.0, 0.0]]) |
---|
[6153] | 899 | assert num.allclose(results.get_attributes(attribute_name='elevation'), |
---|
| 900 | [10.0, 0.0, 10.4]) |
---|
| 901 | assert num.allclose(results.get_attributes(attribute_name='speed'), |
---|
| 902 | [0.0, 10.0, 40.0]) |
---|
[6086] | 903 | |
---|
| 904 | # Blocking |
---|
| 905 | geo_list = [] |
---|
| 906 | for i in results: |
---|
[6304] | 907 | geo_list.append(i) |
---|
[6153] | 908 | assert num.allclose(geo_list[0].get_data_points(), |
---|
| 909 | [[1.0, 0.0],[0.0, 1.0]]) |
---|
| 910 | assert num.allclose(geo_list[0].get_attributes(attribute_name='elevation'), |
---|
| 911 | [10.0, 0.0]) |
---|
| 912 | assert num.allclose(geo_list[1].get_data_points(), |
---|
[6304] | 913 | [[1.0, 0.0]]) |
---|
[6153] | 914 | assert num.allclose(geo_list[1].get_attributes(attribute_name='elevation'), |
---|
| 915 | [10.4]) |
---|
[6086] | 916 | |
---|
[6304] | 917 | os.remove(fileName) |
---|
| 918 | os.remove(pts_file) |
---|
| 919 | |
---|
[6086] | 920 | def verbose_test_load_pts_blocking(self): |
---|
[6304] | 921 | fileName = tempfile.mktemp('.txt') |
---|
| 922 | file = open(fileName, 'w') |
---|
| 923 | file.write(' x,y, elevation , speed \n\ |
---|
[6086] | 924 | 1.0, 0.0, 10.0, 0.0\n\ |
---|
| 925 | 0.0, 1.0, 0.0, 10.0\n\ |
---|
| 926 | 1.0, 0.0, 10.0, 0.0\n\ |
---|
| 927 | 0.0, 1.0, 0.0, 10.0\n\ |
---|
| 928 | 1.0, 0.0, 10.0, 0.0\n\ |
---|
| 929 | 0.0, 1.0, 0.0, 10.0\n\ |
---|
| 930 | 1.0, 0.0, 10.0, 0.0\n\ |
---|
| 931 | 0.0, 1.0, 0.0, 10.0\n\ |
---|
| 932 | 1.0, 0.0, 10.0, 0.0\n\ |
---|
| 933 | 0.0, 1.0, 0.0, 10.0\n\ |
---|
| 934 | 1.0, 0.0, 10.0, 0.0\n\ |
---|
| 935 | 0.0, 1.0, 0.0, 10.0\n\ |
---|
| 936 | 1.0, 0.0, 10.0, 0.0\n\ |
---|
| 937 | 0.0, 1.0, 0.0, 10.0\n\ |
---|
| 938 | 1.0, 0.0, 10.0, 0.0\n\ |
---|
| 939 | 0.0, 1.0, 0.0, 10.0\n\ |
---|
| 940 | 1.0, 0.0, 10.0, 0.0\n\ |
---|
| 941 | 0.0, 1.0, 0.0, 10.0\n\ |
---|
| 942 | 1.0, 0.0, 10.0, 0.0\n\ |
---|
| 943 | 0.0, 1.0, 0.0, 10.0\n\ |
---|
| 944 | 1.0, 0.0, 10.0, 0.0\n\ |
---|
| 945 | 0.0, 1.0, 0.0, 10.0\n\ |
---|
| 946 | 1.0, 0.0, 10.0, 0.0\n\ |
---|
| 947 | 0.0, 1.0, 0.0, 10.0\n\ |
---|
| 948 | 1.0, 0.0, 10.0, 0.0\n\ |
---|
| 949 | 0.0, 1.0, 0.0, 10.0\n\ |
---|
| 950 | 1.0, 0.0, 10.0, 0.0\n\ |
---|
| 951 | 0.0, 1.0, 0.0, 10.0\n\ |
---|
| 952 | 1.0, 0.0, 10.0, 0.0\n\ |
---|
| 953 | 0.0, 1.0, 0.0, 10.0\n\ |
---|
| 954 | 1.0, 0.0, 10.0, 0.0\n\ |
---|
| 955 | 0.0, 1.0, 0.0, 10.0\n\ |
---|
| 956 | 1.0, 0.0, 10.0, 0.0\n\ |
---|
| 957 | 0.0, 1.0, 0.0, 10.0\n\ |
---|
| 958 | 1.0, 0.0, 10.0, 0.0\n\ |
---|
| 959 | 0.0, 1.0, 0.0, 10.0\n\ |
---|
| 960 | 1.0, 0.0, 10.0, 0.0\n\ |
---|
| 961 | 0.0, 1.0, 0.0, 10.0\n\ |
---|
| 962 | 1.0, 0.0, 10.0, 0.0\n\ |
---|
| 963 | 0.0, 1.0, 0.0, 10.0\n\ |
---|
| 964 | 1.0, 0.0, 10.0, 0.0\n\ |
---|
| 965 | 0.0, 1.0, 0.0, 10.0\n\ |
---|
| 966 | 1.0, 0.0, 10.0, 0.0\n\ |
---|
| 967 | 0.0, 1.0, 0.0, 10.0\n\ |
---|
| 968 | 1.0, 0.0, 10.0, 0.0\n\ |
---|
| 969 | 0.0, 1.0, 0.0, 10.0\n\ |
---|
[6304] | 970 | 1.0, 0.0 ,10.4, 40.0\n') |
---|
[6086] | 971 | file.close() |
---|
| 972 | |
---|
[6304] | 973 | pts_file = tempfile.mktemp('.pts') |
---|
[6086] | 974 | convert = Geospatial_data(fileName) |
---|
| 975 | convert.export_points_file(pts_file) |
---|
| 976 | results = Geospatial_data(pts_file, max_read_lines=2, verbose=True) |
---|
| 977 | |
---|
| 978 | # Blocking |
---|
| 979 | geo_list = [] |
---|
| 980 | for i in results: |
---|
[6304] | 981 | geo_list.append(i) |
---|
[6153] | 982 | assert num.allclose(geo_list[0].get_data_points(), |
---|
[6304] | 983 | [[1.0, 0.0], [0.0, 1.0]]) |
---|
[6153] | 984 | assert num.allclose(geo_list[0].get_attributes(attribute_name='elevation'), |
---|
| 985 | [10.0, 0.0]) |
---|
| 986 | assert num.allclose(geo_list[1].get_data_points(), |
---|
[6304] | 987 | [[1.0, 0.0],[0.0, 1.0] ]) |
---|
[6153] | 988 | assert num.allclose(geo_list[1].get_attributes(attribute_name='elevation'), |
---|
| 989 | [10.0, 0.0]) |
---|
[6304] | 990 | |
---|
| 991 | os.remove(fileName) |
---|
[6086] | 992 | os.remove(pts_file) |
---|
| 993 | |
---|
| 994 | def test_new_export_pts_file(self): |
---|
| 995 | att_dict = {} |
---|
[6153] | 996 | pointlist = num.array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]]) |
---|
| 997 | att_dict['elevation'] = num.array([10.1, 0.0, 10.4]) |
---|
| 998 | att_dict['brightness'] = num.array([10.0, 1.0, 10.4]) |
---|
[6304] | 999 | |
---|
| 1000 | fileName = tempfile.mktemp('.pts') |
---|
[6086] | 1001 | G = Geospatial_data(pointlist, att_dict) |
---|
| 1002 | G.export_points_file(fileName) |
---|
| 1003 | results = Geospatial_data(file_name = fileName) |
---|
[6304] | 1004 | os.remove(fileName) |
---|
[6086] | 1005 | |
---|
[6304] | 1006 | assert num.allclose(results.get_data_points(), |
---|
| 1007 | [[1.0, 0.0], [0.0, 1.0], [1.0, 0.0]]) |
---|
| 1008 | assert num.allclose(results.get_attributes(attribute_name='elevation'), |
---|
| 1009 | [10.1, 0.0, 10.4]) |
---|
[6086] | 1010 | answer = [10.0, 1.0, 10.4] |
---|
[6304] | 1011 | assert num.allclose(results.get_attributes(attribute_name='brightness'), |
---|
| 1012 | answer) |
---|
[6086] | 1013 | |
---|
| 1014 | def test_new_export_absolute_pts_file(self): |
---|
| 1015 | att_dict = {} |
---|
[6304] | 1016 | pointlist = num.array([[1.0, 0.0], [0.0, 1.0], [1.0, 0.0]]) |
---|
[6153] | 1017 | att_dict['elevation'] = num.array([10.1, 0.0, 10.4]) |
---|
| 1018 | att_dict['brightness'] = num.array([10.0, 1.0, 10.4]) |
---|
[6086] | 1019 | geo_ref = Geo_reference(50, 25, 55) |
---|
[6304] | 1020 | |
---|
| 1021 | fileName = tempfile.mktemp('.pts') |
---|
[6086] | 1022 | G = Geospatial_data(pointlist, att_dict, geo_ref) |
---|
| 1023 | G.export_points_file(fileName, absolute=True) |
---|
| 1024 | results = Geospatial_data(file_name = fileName) |
---|
[6304] | 1025 | os.remove(fileName) |
---|
[6086] | 1026 | |
---|
[6410] | 1027 | msg = ('results.get_data_points()=\n%s\nbut G.get_data_points(True)=\n%s' |
---|
[6304] | 1028 | % (str(results.get_data_points()), str(G.get_data_points(True)))) |
---|
| 1029 | assert num.allclose(results.get_data_points(), |
---|
| 1030 | G.get_data_points(True)), msg |
---|
| 1031 | msg = ("results.get_attributes(attribute_name='elevation')=%s" |
---|
| 1032 | % str(results.get_attributes(attribute_name='elevation'))) |
---|
| 1033 | assert num.allclose(results.get_attributes(attribute_name='elevation'), |
---|
| 1034 | [10.1, 0.0, 10.4]), msg |
---|
[6086] | 1035 | answer = [10.0, 1.0, 10.4] |
---|
[6304] | 1036 | msg = ("results.get_attributes(attribute_name='brightness')=%s, " |
---|
| 1037 | 'answer=%s' % |
---|
| 1038 | (str(results.get_attributes(attribute_name='brightness')), |
---|
| 1039 | str(answer))) |
---|
| 1040 | assert num.allclose(results.get_attributes(attribute_name='brightness'), |
---|
| 1041 | answer), msg |
---|
[6086] | 1042 | |
---|
| 1043 | def test_loadpts(self): |
---|
| 1044 | from Scientific.IO.NetCDF import NetCDFFile |
---|
| 1045 | |
---|
[6304] | 1046 | fileName = tempfile.mktemp('.pts') |
---|
[6086] | 1047 | # NetCDF file definition |
---|
| 1048 | outfile = NetCDFFile(fileName, netcdf_mode_w) |
---|
[6304] | 1049 | |
---|
[6086] | 1050 | # dimension definitions |
---|
[6304] | 1051 | outfile.createDimension('number_of_points', 3) |
---|
| 1052 | outfile.createDimension('number_of_dimensions', 2) # This is 2d data |
---|
| 1053 | |
---|
[6086] | 1054 | # variable definitions |
---|
[6304] | 1055 | outfile.createVariable('points', netcdf_float, ('number_of_points', |
---|
| 1056 | 'number_of_dimensions')) |
---|
| 1057 | outfile.createVariable('elevation', netcdf_float, ('number_of_points',)) |
---|
| 1058 | |
---|
[6086] | 1059 | # Get handles to the variables |
---|
| 1060 | points = outfile.variables['points'] |
---|
| 1061 | elevation = outfile.variables['elevation'] |
---|
[6304] | 1062 | |
---|
[6086] | 1063 | points[0, :] = [1.0,0.0] |
---|
[6304] | 1064 | elevation[0] = 10.0 |
---|
[6086] | 1065 | points[1, :] = [0.0,1.0] |
---|
[6304] | 1066 | elevation[1] = 0.0 |
---|
[6086] | 1067 | points[2, :] = [1.0,0.0] |
---|
[6304] | 1068 | elevation[2] = 10.4 |
---|
[6086] | 1069 | |
---|
| 1070 | outfile.close() |
---|
[6304] | 1071 | |
---|
[6086] | 1072 | results = Geospatial_data(file_name = fileName) |
---|
| 1073 | os.remove(fileName) |
---|
[6304] | 1074 | |
---|
| 1075 | answer = [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]] |
---|
| 1076 | assert num.allclose(results.get_data_points(), |
---|
| 1077 | [[1.0, 0.0], [0.0, 1.0], [1.0, 0.0]]) |
---|
| 1078 | assert num.allclose(results.get_attributes(attribute_name='elevation'), |
---|
| 1079 | [10.0, 0.0, 10.4]) |
---|
| 1080 | |
---|
[6086] | 1081 | def test_writepts(self): |
---|
[6304] | 1082 | '''Test that storage of x,y,attributes works''' |
---|
| 1083 | |
---|
[6086] | 1084 | att_dict = {} |
---|
[6304] | 1085 | pointlist = num.array([[1.0, 0.0], [0.0, 1.0], [1.0, 0.0]]) |
---|
[6153] | 1086 | att_dict['elevation'] = num.array([10.0, 0.0, 10.4]) |
---|
| 1087 | att_dict['brightness'] = num.array([10.0, 0.0, 10.4]) |
---|
[6304] | 1088 | geo_reference=Geo_reference(56, 1.9, 1.9) |
---|
[6086] | 1089 | |
---|
| 1090 | # Test pts format |
---|
[6304] | 1091 | fileName = tempfile.mktemp('.pts') |
---|
[6086] | 1092 | G = Geospatial_data(pointlist, att_dict, geo_reference) |
---|
| 1093 | G.export_points_file(fileName, False) |
---|
| 1094 | results = Geospatial_data(file_name=fileName) |
---|
| 1095 | os.remove(fileName) |
---|
| 1096 | |
---|
[6304] | 1097 | assert num.allclose(results.get_data_points(False), |
---|
| 1098 | [[1.0, 0.0], [0.0, 1.0], [1.0, 0.0]]) |
---|
| 1099 | assert num.allclose(results.get_attributes('elevation'), |
---|
| 1100 | [10.0, 0.0, 10.4]) |
---|
[6086] | 1101 | answer = [10.0, 0.0, 10.4] |
---|
[6153] | 1102 | assert num.allclose(results.get_attributes('brightness'), answer) |
---|
[6086] | 1103 | self.failUnless(geo_reference == geo_reference, |
---|
[6304] | 1104 | 'test_writepts failed. Test geo_reference') |
---|
[6086] | 1105 | |
---|
| 1106 | def test_write_csv_attributes(self): |
---|
[6304] | 1107 | '''Test that storage of x,y,attributes works''' |
---|
| 1108 | |
---|
[6086] | 1109 | att_dict = {} |
---|
[6304] | 1110 | pointlist = num.array([[1.0, 0.0], [0.0, 1.0], [1.0, 0.0]]) |
---|
[6153] | 1111 | att_dict['elevation'] = num.array([10.0, 0.0, 10.4]) |
---|
| 1112 | att_dict['brightness'] = num.array([10.0, 0.0, 10.4]) |
---|
[6304] | 1113 | geo_reference=Geo_reference(56, 0, 0) |
---|
| 1114 | |
---|
[6086] | 1115 | # Test txt format |
---|
[6304] | 1116 | fileName = tempfile.mktemp('.txt') |
---|
[6086] | 1117 | G = Geospatial_data(pointlist, att_dict, geo_reference) |
---|
| 1118 | G.export_points_file(fileName) |
---|
| 1119 | results = Geospatial_data(file_name=fileName) |
---|
| 1120 | os.remove(fileName) |
---|
[6304] | 1121 | |
---|
| 1122 | assert num.allclose(results.get_data_points(False), |
---|
| 1123 | [[1.0, 0.0], [0.0, 1.0], [1.0, 0.0]]) |
---|
| 1124 | assert num.allclose(results.get_attributes('elevation'), |
---|
| 1125 | [10.0, 0.0, 10.4]) |
---|
[6086] | 1126 | answer = [10.0, 0.0, 10.4] |
---|
[6153] | 1127 | assert num.allclose(results.get_attributes('brightness'), answer) |
---|
[6304] | 1128 | |
---|
[6086] | 1129 | def test_write_csv_attributes_lat_long(self): |
---|
[6304] | 1130 | '''Test that storage of x,y,attributes works''' |
---|
| 1131 | |
---|
[6086] | 1132 | att_dict = {} |
---|
[6304] | 1133 | pointlist = num.array([[-21.5, 114.5], [-21.6, 114.5], [-21.7, 114.5]]) |
---|
[6153] | 1134 | att_dict['elevation'] = num.array([10.0, 0.0, 10.4]) |
---|
| 1135 | att_dict['brightness'] = num.array([10.0, 0.0, 10.4]) |
---|
[6304] | 1136 | |
---|
[6086] | 1137 | # Test txt format |
---|
[6304] | 1138 | fileName = tempfile.mktemp('.txt') |
---|
[6086] | 1139 | G = Geospatial_data(pointlist, att_dict, points_are_lats_longs=True) |
---|
| 1140 | G.export_points_file(fileName, as_lat_long=True) |
---|
| 1141 | results = Geospatial_data(file_name=fileName) |
---|
| 1142 | os.remove(fileName) |
---|
[6304] | 1143 | |
---|
[6153] | 1144 | assert num.allclose(results.get_data_points(False, as_lat_long=True), |
---|
[6304] | 1145 | pointlist) |
---|
| 1146 | assert num.allclose(results.get_attributes('elevation'), |
---|
| 1147 | [10.0, 0.0, 10.4]) |
---|
[6086] | 1148 | answer = [10.0, 0.0, 10.4] |
---|
[6153] | 1149 | assert num.allclose(results.get_attributes('brightness'), answer) |
---|
[6304] | 1150 | |
---|
[6086] | 1151 | def test_writepts_no_attributes(self): |
---|
[6304] | 1152 | '''Test that storage of x,y alone works''' |
---|
[6086] | 1153 | |
---|
| 1154 | att_dict = {} |
---|
[6304] | 1155 | pointlist = num.array([[1.0, 0.0], [0.0, 1.0], [1.0, 0.0]]) |
---|
| 1156 | geo_reference=Geo_reference(56, 1.9, 1.9) |
---|
[6086] | 1157 | |
---|
| 1158 | # Test pts format |
---|
[6304] | 1159 | fileName = tempfile.mktemp('.pts') |
---|
[6086] | 1160 | G = Geospatial_data(pointlist, None, geo_reference) |
---|
| 1161 | G.export_points_file(fileName, False) |
---|
| 1162 | results = Geospatial_data(file_name=fileName) |
---|
| 1163 | os.remove(fileName) |
---|
| 1164 | |
---|
[6304] | 1165 | assert num.allclose(results.get_data_points(False), |
---|
| 1166 | [[1.0, 0.0], [0.0, 1.0], [1.0, 0.0]]) |
---|
[6086] | 1167 | self.failUnless(geo_reference == geo_reference, |
---|
[6304] | 1168 | 'test_writepts failed. Test geo_reference') |
---|
| 1169 | |
---|
[6086] | 1170 | def test_write_csv_no_attributes(self): |
---|
[6304] | 1171 | '''Test that storage of x,y alone works''' |
---|
| 1172 | |
---|
[6086] | 1173 | att_dict = {} |
---|
[6304] | 1174 | pointlist = num.array([[1.0, 0.0], [0.0, 1.0], [1.0, 0.0]]) |
---|
[6086] | 1175 | geo_reference=Geo_reference(56,0,0) |
---|
[6304] | 1176 | |
---|
[6086] | 1177 | # Test format |
---|
[6304] | 1178 | fileName = tempfile.mktemp('.txt') |
---|
[6086] | 1179 | G = Geospatial_data(pointlist, None, geo_reference) |
---|
| 1180 | G.export_points_file(fileName) |
---|
| 1181 | results = Geospatial_data(file_name=fileName) |
---|
| 1182 | os.remove(fileName) |
---|
| 1183 | |
---|
[6304] | 1184 | assert num.allclose(results.get_data_points(False), |
---|
| 1185 | [[1.0, 0.0], [0.0, 1.0], [1.0, 0.0]]) |
---|
[6086] | 1186 | |
---|
[6304] | 1187 | ################################################################################ |
---|
| 1188 | # Check bad PTS files. |
---|
| 1189 | ################################################################################ |
---|
| 1190 | |
---|
[6086] | 1191 | def test_load_bad_no_file_pts(self): |
---|
[6304] | 1192 | fileName = tempfile.mktemp('.pts') |
---|
[6086] | 1193 | try: |
---|
[6304] | 1194 | results = Geospatial_data(file_name=fileName) |
---|
[6086] | 1195 | except IOError: |
---|
| 1196 | pass |
---|
| 1197 | else: |
---|
| 1198 | msg = 'imaginary file did not raise error!' |
---|
[6304] | 1199 | raise Exception, msg |
---|
[6086] | 1200 | |
---|
| 1201 | def test_create_from_pts_file(self): |
---|
| 1202 | from Scientific.IO.NetCDF import NetCDFFile |
---|
| 1203 | |
---|
[6304] | 1204 | # NetCDF file definition |
---|
[6086] | 1205 | FN = 'test_points.pts' |
---|
| 1206 | outfile = NetCDFFile(FN, netcdf_mode_w) |
---|
[6304] | 1207 | |
---|
[6086] | 1208 | # dimension definitions |
---|
[6304] | 1209 | outfile.createDimension('number_of_points', 3) |
---|
| 1210 | outfile.createDimension('number_of_dimensions', 2) # This is 2d data |
---|
| 1211 | |
---|
[6086] | 1212 | # variable definitions |
---|
[6304] | 1213 | outfile.createVariable('points', netcdf_float, ('number_of_points', |
---|
| 1214 | 'number_of_dimensions')) |
---|
| 1215 | outfile.createVariable('elevation', netcdf_float, ('number_of_points',)) |
---|
| 1216 | |
---|
[6086] | 1217 | # Get handles to the variables |
---|
| 1218 | points = outfile.variables['points'] |
---|
| 1219 | elevation = outfile.variables['elevation'] |
---|
[6304] | 1220 | |
---|
[6086] | 1221 | points[0, :] = [1.0,0.0] |
---|
[6304] | 1222 | elevation[0] = 10.0 |
---|
[6086] | 1223 | points[1, :] = [0.0,1.0] |
---|
[6304] | 1224 | elevation[1] = 0.0 |
---|
[6086] | 1225 | points[2, :] = [1.0,0.0] |
---|
[6304] | 1226 | elevation[2] = 10.4 |
---|
[6086] | 1227 | |
---|
| 1228 | outfile.close() |
---|
| 1229 | |
---|
| 1230 | G = Geospatial_data(file_name = FN) |
---|
| 1231 | |
---|
[6153] | 1232 | assert num.allclose(G.get_geo_reference().get_xllcorner(), 0.0) |
---|
| 1233 | assert num.allclose(G.get_geo_reference().get_yllcorner(), 0.0) |
---|
[6304] | 1234 | assert num.allclose(G.get_data_points(), |
---|
| 1235 | [[1.0, 0.0], [0.0, 1.0], [1.0, 0.0]]) |
---|
[6153] | 1236 | assert num.allclose(G.get_attributes(), [10.0, 0.0, 10.4]) |
---|
[6086] | 1237 | os.remove(FN) |
---|
| 1238 | |
---|
| 1239 | def test_create_from_pts_file_with_geo(self): |
---|
[6304] | 1240 | '''Test if Geospatial data is correctly instantiated from a pts file.''' |
---|
| 1241 | |
---|
[6086] | 1242 | from Scientific.IO.NetCDF import NetCDFFile |
---|
| 1243 | |
---|
[6304] | 1244 | # NetCDF file definition |
---|
[6086] | 1245 | FN = 'test_points.pts' |
---|
| 1246 | outfile = NetCDFFile(FN, netcdf_mode_w) |
---|
| 1247 | |
---|
| 1248 | # Make up an arbitrary georef |
---|
| 1249 | xll = 0.1 |
---|
| 1250 | yll = 20 |
---|
| 1251 | geo_reference=Geo_reference(56, xll, yll) |
---|
| 1252 | geo_reference.write_NetCDF(outfile) |
---|
| 1253 | |
---|
| 1254 | # dimension definitions |
---|
[6304] | 1255 | outfile.createDimension('number_of_points', 3) |
---|
| 1256 | outfile.createDimension('number_of_dimensions', 2) # This is 2d data |
---|
| 1257 | |
---|
[6086] | 1258 | # variable definitions |
---|
[6304] | 1259 | outfile.createVariable('points', netcdf_float, ('number_of_points', |
---|
| 1260 | 'number_of_dimensions')) |
---|
| 1261 | outfile.createVariable('elevation', netcdf_float, ('number_of_points',)) |
---|
| 1262 | |
---|
[6086] | 1263 | # Get handles to the variables |
---|
| 1264 | points = outfile.variables['points'] |
---|
| 1265 | elevation = outfile.variables['elevation'] |
---|
| 1266 | |
---|
| 1267 | points[0, :] = [1.0,0.0] |
---|
[6304] | 1268 | elevation[0] = 10.0 |
---|
[6086] | 1269 | points[1, :] = [0.0,1.0] |
---|
[6304] | 1270 | elevation[1] = 0.0 |
---|
[6086] | 1271 | points[2, :] = [1.0,0.0] |
---|
[6304] | 1272 | elevation[2] = 10.4 |
---|
[6086] | 1273 | |
---|
| 1274 | outfile.close() |
---|
| 1275 | |
---|
| 1276 | G = Geospatial_data(file_name = FN) |
---|
| 1277 | |
---|
[6153] | 1278 | assert num.allclose(G.get_geo_reference().get_xllcorner(), xll) |
---|
| 1279 | assert num.allclose(G.get_geo_reference().get_yllcorner(), yll) |
---|
| 1280 | assert num.allclose(G.get_data_points(), [[1.0+xll, 0.0+yll], |
---|
| 1281 | [0.0+xll, 1.0+yll], |
---|
| 1282 | [1.0+xll, 0.0+yll]]) |
---|
| 1283 | assert num.allclose(G.get_attributes(), [10.0, 0.0, 10.4]) |
---|
[6304] | 1284 | |
---|
[6086] | 1285 | os.remove(FN) |
---|
| 1286 | |
---|
| 1287 | def test_add_(self): |
---|
| 1288 | '''test_add_(self): |
---|
| 1289 | adds an txt and pts files, reads the files and adds them |
---|
[6304] | 1290 | checking results are correct |
---|
[6086] | 1291 | ''' |
---|
[6304] | 1292 | |
---|
[6086] | 1293 | # create files |
---|
| 1294 | att_dict1 = {} |
---|
[6304] | 1295 | pointlist1 = num.array([[1.0, 0.0], [0.0, 1.0], [1.0, 0.0]]) |
---|
[6153] | 1296 | att_dict1['elevation'] = num.array([-10.0, 0.0, 10.4]) |
---|
| 1297 | att_dict1['brightness'] = num.array([10.0, 0.0, 10.4]) |
---|
[6086] | 1298 | geo_reference1 = Geo_reference(56, 2.0, 1.0) |
---|
[6304] | 1299 | |
---|
[6086] | 1300 | att_dict2 = {} |
---|
[6304] | 1301 | pointlist2 = num.array([[2.0, 1.0], [1.0, 2.0], [2.0, 1.0]]) |
---|
[6153] | 1302 | att_dict2['elevation'] = num.array([1.0, 15.0, 1.4]) |
---|
| 1303 | att_dict2['brightness'] = num.array([14.0, 1.0, -12.4]) |
---|
[6304] | 1304 | geo_reference2 = Geo_reference(56, 1.0, 2.0) |
---|
[6086] | 1305 | |
---|
| 1306 | G1 = Geospatial_data(pointlist1, att_dict1, geo_reference1) |
---|
| 1307 | G2 = Geospatial_data(pointlist2, att_dict2, geo_reference2) |
---|
| 1308 | |
---|
[6304] | 1309 | fileName1 = tempfile.mktemp('.txt') |
---|
| 1310 | fileName2 = tempfile.mktemp('.pts') |
---|
| 1311 | |
---|
| 1312 | # makes files |
---|
[6086] | 1313 | G1.export_points_file(fileName1) |
---|
| 1314 | G2.export_points_file(fileName2) |
---|
[6304] | 1315 | |
---|
[6086] | 1316 | # add files |
---|
[6304] | 1317 | G3 = Geospatial_data(file_name=fileName1) |
---|
| 1318 | G4 = Geospatial_data(file_name=fileName2) |
---|
[6086] | 1319 | G = G3 + G4 |
---|
| 1320 | |
---|
| 1321 | #read results |
---|
[6153] | 1322 | assert num.allclose(G.get_data_points(), |
---|
| 1323 | [[ 3.0, 1.0], [ 2.0, 2.0], |
---|
| 1324 | [ 3.0, 1.0], [ 3.0, 3.0], |
---|
| 1325 | [ 2.0, 4.0], [ 3.0, 3.0]]) |
---|
| 1326 | assert num.allclose(G.get_attributes(attribute_name='elevation'), |
---|
| 1327 | [-10.0, 0.0, 10.4, 1.0, 15.0, 1.4]) |
---|
[6086] | 1328 | answer = [10.0, 0.0, 10.4, 14.0, 1.0, -12.4] |
---|
[6304] | 1329 | assert num.allclose(G.get_attributes(attribute_name='brightness'), |
---|
| 1330 | answer) |
---|
[6086] | 1331 | self.failUnless(G.get_geo_reference() == geo_reference1, |
---|
[6304] | 1332 | 'test_writepts failed. Test geo_reference') |
---|
| 1333 | |
---|
[6086] | 1334 | os.remove(fileName1) |
---|
| 1335 | os.remove(fileName2) |
---|
[6304] | 1336 | |
---|
[6086] | 1337 | def test_ensure_absolute(self): |
---|
[6304] | 1338 | points = [[2.0, 0.0], [1.0, 1.0], |
---|
| 1339 | [2.0, 0.0], [2.0, 2.0], |
---|
| 1340 | [1.0, 3.0], [2.0, 2.0]] |
---|
[6086] | 1341 | new_points = ensure_absolute(points) |
---|
[6304] | 1342 | |
---|
[6153] | 1343 | assert num.allclose(new_points, points) |
---|
[6086] | 1344 | |
---|
[6304] | 1345 | points = num.array([[2.0, 0.0], [1.0, 1.0], |
---|
| 1346 | [2.0, 0.0], [2.0, 2.0], |
---|
| 1347 | [1.0, 3.0], [2.0, 2.0]]) |
---|
[6086] | 1348 | new_points = ensure_absolute(points) |
---|
[6304] | 1349 | |
---|
[6153] | 1350 | assert num.allclose(new_points, points) |
---|
[6086] | 1351 | |
---|
[6304] | 1352 | ab_points = num.array([[2.0, 0.0], [1.0, 1.0], |
---|
| 1353 | [2.0, 0.0], [2.0, 2.0], |
---|
| 1354 | [1.0, 3.0], [2.0, 2.0]]) |
---|
| 1355 | |
---|
| 1356 | mesh_origin = (56, 290000, 618000) # zone, easting, northing |
---|
| 1357 | data_points = num.zeros((ab_points.shape), num.float) |
---|
| 1358 | |
---|
[6086] | 1359 | #Shift datapoints according to new origins |
---|
| 1360 | for k in range(len(ab_points)): |
---|
| 1361 | data_points[k][0] = ab_points[k][0] - mesh_origin[1] |
---|
| 1362 | data_points[k][1] = ab_points[k][1] - mesh_origin[2] |
---|
[6304] | 1363 | new_points = ensure_absolute(data_points, geo_reference=mesh_origin) |
---|
| 1364 | |
---|
[6153] | 1365 | assert num.allclose(new_points, ab_points) |
---|
[6086] | 1366 | |
---|
| 1367 | geo = Geo_reference(56,67,-56) |
---|
[6304] | 1368 | data_points = geo.change_points_geo_ref(ab_points) |
---|
| 1369 | new_points = ensure_absolute(data_points, geo_reference=geo) |
---|
[6086] | 1370 | |
---|
[6153] | 1371 | assert num.allclose(new_points, ab_points) |
---|
[6086] | 1372 | |
---|
| 1373 | geo_reference = Geo_reference(56, 100, 200) |
---|
| 1374 | ab_points = [[1.0, 2.1], [3.0, 5.3]] |
---|
| 1375 | points = geo_reference.change_points_geo_ref(ab_points) |
---|
| 1376 | attributes = [2, 4] |
---|
[6304] | 1377 | G = Geospatial_data(points, attributes, geo_reference=geo_reference) |
---|
[6086] | 1378 | new_points = ensure_absolute(G) |
---|
[6304] | 1379 | |
---|
[6153] | 1380 | assert num.allclose(new_points, ab_points) |
---|
[6086] | 1381 | |
---|
| 1382 | def test_ensure_geospatial(self): |
---|
[6304] | 1383 | points = [[2.0, 0.0], [1.0, 1.0], |
---|
| 1384 | [2.0, 0.0], [2.0, 2.0], |
---|
| 1385 | [1.0, 3.0], [2.0, 2.0]] |
---|
[6086] | 1386 | new_points = ensure_geospatial(points) |
---|
[6304] | 1387 | assert num.allclose(new_points.get_data_points(absolute=True), points) |
---|
[6086] | 1388 | |
---|
[6304] | 1389 | points = num.array([[2.0, 0.0], [1.0, 1.0], |
---|
| 1390 | [2.0, 0.0], [2.0, 2.0], |
---|
| 1391 | [1.0, 3.0], [2.0, 2.0]]) |
---|
[6086] | 1392 | new_points = ensure_geospatial(points) |
---|
[6304] | 1393 | assert num.allclose(new_points.get_data_points(absolute=True), points) |
---|
| 1394 | |
---|
[6153] | 1395 | ab_points = num.array([[2.0, 0.0],[1.0, 1.0], |
---|
| 1396 | [2.0, 0.0],[2.0, 2.0], |
---|
| 1397 | [1.0, 3.0],[2.0, 2.0]]) |
---|
[6304] | 1398 | mesh_origin = (56, 290000, 618000) # zone, easting, northing |
---|
| 1399 | data_points = num.zeros((ab_points.shape), num.float) |
---|
[6086] | 1400 | |
---|
| 1401 | #Shift datapoints according to new origins |
---|
| 1402 | for k in range(len(ab_points)): |
---|
| 1403 | data_points[k][0] = ab_points[k][0] - mesh_origin[1] |
---|
| 1404 | data_points[k][1] = ab_points[k][1] - mesh_origin[2] |
---|
| 1405 | new_geospatial = ensure_geospatial(data_points, |
---|
[6304] | 1406 | geo_reference=mesh_origin) |
---|
[6086] | 1407 | new_points = new_geospatial.get_data_points(absolute=True) |
---|
[6153] | 1408 | assert num.allclose(new_points, ab_points) |
---|
[6086] | 1409 | |
---|
[6304] | 1410 | geo = Geo_reference(56, 67, -56) |
---|
| 1411 | data_points = geo.change_points_geo_ref(ab_points) |
---|
| 1412 | new_geospatial = ensure_geospatial(data_points, geo_reference=geo) |
---|
[6086] | 1413 | new_points = new_geospatial.get_data_points(absolute=True) |
---|
[6441] | 1414 | msg = ('new_points=\n%s\nab_points=\n%s' |
---|
| 1415 | % (str(new_points), str(ab_points))) |
---|
| 1416 | assert num.allclose(new_points, ab_points), msg |
---|
[6086] | 1417 | |
---|
| 1418 | geo_reference = Geo_reference(56, 100, 200) |
---|
| 1419 | ab_points = [[1.0, 2.1], [3.0, 5.3]] |
---|
| 1420 | points = geo_reference.change_points_geo_ref(ab_points) |
---|
| 1421 | attributes = [2, 4] |
---|
[6304] | 1422 | G = Geospatial_data(points, attributes, geo_reference=geo_reference) |
---|
| 1423 | new_geospatial = ensure_geospatial(G) |
---|
[6086] | 1424 | new_points = new_geospatial.get_data_points(absolute=True) |
---|
[6153] | 1425 | assert num.allclose(new_points, ab_points) |
---|
[6304] | 1426 | |
---|
[6086] | 1427 | def test_isinstance(self): |
---|
[6304] | 1428 | fileName = tempfile.mktemp('.csv') |
---|
| 1429 | file = open(fileName, 'w') |
---|
| 1430 | file.write('x,y, elevation , speed \n\ |
---|
[6086] | 1431 | 1.0, 0.0, 10.0, 0.0\n\ |
---|
| 1432 | 0.0, 1.0, 0.0, 10.0\n\ |
---|
[6304] | 1433 | 1.0, 0.0, 10.4, 40.0\n') |
---|
[6086] | 1434 | file.close() |
---|
| 1435 | |
---|
| 1436 | results = Geospatial_data(fileName) |
---|
[6153] | 1437 | assert num.allclose(results.get_data_points(absolute=True), |
---|
[6304] | 1438 | [[1.0, 0.0], [0.0, 1.0], [1.0, 0.0]]) |
---|
[6153] | 1439 | assert num.allclose(results.get_attributes(attribute_name='elevation'), |
---|
| 1440 | [10.0, 0.0, 10.4]) |
---|
| 1441 | assert num.allclose(results.get_attributes(attribute_name='speed'), |
---|
| 1442 | [0.0, 10.0, 40.0]) |
---|
[6086] | 1443 | |
---|
| 1444 | os.remove(fileName) |
---|
| 1445 | |
---|
| 1446 | def test_no_constructors(self): |
---|
| 1447 | try: |
---|
| 1448 | G = Geospatial_data() |
---|
| 1449 | except ValueError: |
---|
| 1450 | pass |
---|
| 1451 | else: |
---|
| 1452 | msg = 'Instance must have a filename or data points' |
---|
[6304] | 1453 | raise Exception, msg |
---|
[6086] | 1454 | |
---|
| 1455 | def test_load_csv_lat_long(self): |
---|
[6304] | 1456 | '''comma delimited''' |
---|
[6086] | 1457 | |
---|
[6304] | 1458 | fileName = tempfile.mktemp('.csv') |
---|
| 1459 | file = open(fileName, 'w') |
---|
| 1460 | file.write('long,lat, elevation, yeah \n\ |
---|
[6086] | 1461 | 150.916666667,-34.50,452.688000, 10\n\ |
---|
[6304] | 1462 | 150.0,-34,459.126000, 10\n') |
---|
[6086] | 1463 | file.close() |
---|
[6304] | 1464 | |
---|
[6086] | 1465 | results = Geospatial_data(fileName) |
---|
| 1466 | os.remove(fileName) |
---|
| 1467 | points = results.get_data_points() |
---|
[6304] | 1468 | |
---|
[6153] | 1469 | assert num.allclose(points[0][0], 308728.009) |
---|
| 1470 | assert num.allclose(points[0][1], 6180432.601) |
---|
| 1471 | assert num.allclose(points[1][0], 222908.705) |
---|
| 1472 | assert num.allclose(points[1][1], 6233785.284) |
---|
[6304] | 1473 | |
---|
[6086] | 1474 | def test_load_csv_lat_longII(self): |
---|
[6304] | 1475 | '''comma delimited''' |
---|
[6086] | 1476 | |
---|
[6304] | 1477 | fileName = tempfile.mktemp('.csv') |
---|
| 1478 | file = open(fileName, 'w') |
---|
| 1479 | file.write('Lati,LONG,z \n\ |
---|
[6086] | 1480 | -34.50,150.916666667,452.688000\n\ |
---|
[6304] | 1481 | -34,150.0,459.126000\n') |
---|
[6086] | 1482 | file.close() |
---|
[6304] | 1483 | |
---|
[6086] | 1484 | results = Geospatial_data(fileName) |
---|
| 1485 | os.remove(fileName) |
---|
| 1486 | points = results.get_data_points() |
---|
[6304] | 1487 | |
---|
[6153] | 1488 | assert num.allclose(points[0][0], 308728.009) |
---|
| 1489 | assert num.allclose(points[0][1], 6180432.601) |
---|
[6304] | 1490 | assert num.allclose(points[1][0], 222908.705) |
---|
[6153] | 1491 | assert num.allclose(points[1][1], 6233785.284) |
---|
[6086] | 1492 | |
---|
| 1493 | def test_load_csv_lat_long_bad(self): |
---|
[6304] | 1494 | '''comma delimited''' |
---|
[6086] | 1495 | |
---|
[6304] | 1496 | fileName = tempfile.mktemp('.csv') |
---|
| 1497 | file = open(fileName, 'w') |
---|
| 1498 | file.write('Lati,LONG,z \n\ |
---|
[6086] | 1499 | -25.0,180.0,452.688000\n\ |
---|
[6304] | 1500 | -34,150.0,459.126000\n') |
---|
[6086] | 1501 | file.close() |
---|
[6304] | 1502 | |
---|
[6086] | 1503 | try: |
---|
| 1504 | results = Geospatial_data(fileName) |
---|
| 1505 | except ANUGAError: |
---|
| 1506 | pass |
---|
| 1507 | else: |
---|
| 1508 | msg = 'Different zones in Geo references not caught.' |
---|
[6304] | 1509 | raise Exception, msg |
---|
| 1510 | |
---|
[6086] | 1511 | os.remove(fileName) |
---|
[6304] | 1512 | |
---|
[6086] | 1513 | def test_lat_long(self): |
---|
[6304] | 1514 | lat_gong = degminsec2decimal_degrees(-34, 30, 0.) |
---|
| 1515 | lon_gong = degminsec2decimal_degrees(150, 55, 0.) |
---|
| 1516 | |
---|
| 1517 | lat_2 = degminsec2decimal_degrees(-34, 00, 0.) |
---|
| 1518 | lon_2 = degminsec2decimal_degrees(150, 00, 0.) |
---|
| 1519 | |
---|
[6086] | 1520 | lats = [lat_gong, lat_2] |
---|
| 1521 | longs = [lon_gong, lon_2] |
---|
| 1522 | gsd = Geospatial_data(latitudes=lats, longitudes=longs) |
---|
| 1523 | |
---|
| 1524 | points = gsd.get_data_points(absolute=True) |
---|
[6304] | 1525 | |
---|
[6153] | 1526 | assert num.allclose(points[0][0], 308728.009) |
---|
| 1527 | assert num.allclose(points[0][1], 6180432.601) |
---|
[6304] | 1528 | assert num.allclose(points[1][0], 222908.705) |
---|
[6153] | 1529 | assert num.allclose(points[1][1], 6233785.284) |
---|
[6086] | 1530 | self.failUnless(gsd.get_geo_reference().get_zone() == 56, |
---|
| 1531 | 'Bad zone error!') |
---|
[6304] | 1532 | |
---|
[6086] | 1533 | try: |
---|
| 1534 | results = Geospatial_data(latitudes=lats) |
---|
| 1535 | except ValueError: |
---|
| 1536 | pass |
---|
| 1537 | else: |
---|
[6304] | 1538 | self.fail('Error not thrown error!') |
---|
[6086] | 1539 | try: |
---|
| 1540 | results = Geospatial_data(latitudes=lats) |
---|
| 1541 | except ValueError: |
---|
| 1542 | pass |
---|
| 1543 | else: |
---|
[6304] | 1544 | self.fail('Error not thrown error!') |
---|
[6086] | 1545 | try: |
---|
| 1546 | results = Geospatial_data(longitudes=lats) |
---|
| 1547 | except ValueError: |
---|
| 1548 | pass |
---|
| 1549 | else: |
---|
[6304] | 1550 | self.fail('Error not thrown error!') |
---|
[6086] | 1551 | try: |
---|
| 1552 | results = Geospatial_data(latitudes=lats, longitudes=longs, |
---|
| 1553 | geo_reference="p") |
---|
| 1554 | except ValueError: |
---|
| 1555 | pass |
---|
| 1556 | else: |
---|
[6304] | 1557 | self.fail('Error not thrown error!') |
---|
| 1558 | |
---|
[6086] | 1559 | try: |
---|
| 1560 | results = Geospatial_data(latitudes=lats, longitudes=longs, |
---|
| 1561 | data_points=12) |
---|
| 1562 | except ValueError: |
---|
| 1563 | pass |
---|
| 1564 | else: |
---|
[6304] | 1565 | self.fail('Error not thrown error!') |
---|
[6086] | 1566 | |
---|
| 1567 | def test_lat_long2(self): |
---|
[6304] | 1568 | lat_gong = degminsec2decimal_degrees(-34, 30, 0.) |
---|
| 1569 | lon_gong = degminsec2decimal_degrees(150, 55, 0.) |
---|
| 1570 | |
---|
| 1571 | lat_2 = degminsec2decimal_degrees(-34, 00, 0.) |
---|
| 1572 | lon_2 = degminsec2decimal_degrees(150, 00, 0.) |
---|
| 1573 | |
---|
[6086] | 1574 | points = [[lat_gong, lon_gong], [lat_2, lon_2]] |
---|
| 1575 | gsd = Geospatial_data(data_points=points, points_are_lats_longs=True) |
---|
| 1576 | |
---|
| 1577 | points = gsd.get_data_points(absolute=True) |
---|
[6304] | 1578 | |
---|
[6153] | 1579 | assert num.allclose(points[0][0], 308728.009) |
---|
| 1580 | assert num.allclose(points[0][1], 6180432.601) |
---|
[6304] | 1581 | assert num.allclose(points[1][0], 222908.705) |
---|
[6153] | 1582 | assert num.allclose(points[1][1], 6233785.284) |
---|
[6086] | 1583 | self.failUnless(gsd.get_geo_reference().get_zone() == 56, |
---|
| 1584 | 'Bad zone error!') |
---|
| 1585 | |
---|
| 1586 | try: |
---|
| 1587 | results = Geospatial_data(points_are_lats_longs=True) |
---|
| 1588 | except ValueError: |
---|
| 1589 | pass |
---|
| 1590 | else: |
---|
[6304] | 1591 | self.fail('Error not thrown error!') |
---|
[6086] | 1592 | |
---|
[6304] | 1593 | def test_write_urs_file(self): |
---|
| 1594 | lat_gong = degminsec2decimal_degrees(-34, 00, 0) |
---|
| 1595 | lon_gong = degminsec2decimal_degrees(150, 30, 0.) |
---|
[6086] | 1596 | |
---|
[6304] | 1597 | lat_2 = degminsec2decimal_degrees(-34, 00, 1) |
---|
| 1598 | lon_2 = degminsec2decimal_degrees(150, 00, 0.) |
---|
[6086] | 1599 | p1 = (lat_gong, lon_gong) |
---|
| 1600 | p2 = (lat_2, lon_2) |
---|
| 1601 | points = ImmutableSet([p1, p2, p1]) |
---|
| 1602 | gsd = Geospatial_data(data_points=list(points), |
---|
| 1603 | points_are_lats_longs=True) |
---|
[6304] | 1604 | |
---|
| 1605 | fn = tempfile.mktemp('.urs') |
---|
[6086] | 1606 | gsd.export_points_file(fn) |
---|
| 1607 | handle = open(fn) |
---|
| 1608 | lines = handle.readlines() |
---|
[6304] | 1609 | assert lines[0], '2' |
---|
| 1610 | assert lines[1], '-34.0002778 150.0 0' |
---|
| 1611 | assert lines[2], '-34.0 150.5 1' |
---|
[6086] | 1612 | handle.close() |
---|
| 1613 | os.remove(fn) |
---|
[6304] | 1614 | |
---|
[6086] | 1615 | def test_lat_long_set(self): |
---|
[6304] | 1616 | lat_gong = degminsec2decimal_degrees(-34, 30, 0.) |
---|
| 1617 | lon_gong = degminsec2decimal_degrees(150, 55, 0.) |
---|
| 1618 | |
---|
| 1619 | lat_2 = degminsec2decimal_degrees(-34, 00, 0.) |
---|
| 1620 | lon_2 = degminsec2decimal_degrees(150, 00, 0.) |
---|
[6086] | 1621 | p1 = (lat_gong, lon_gong) |
---|
| 1622 | p2 = (lat_2, lon_2) |
---|
| 1623 | points = ImmutableSet([p1, p2, p1]) |
---|
| 1624 | gsd = Geospatial_data(data_points=list(points), |
---|
| 1625 | points_are_lats_longs=True) |
---|
| 1626 | |
---|
| 1627 | points = gsd.get_data_points(absolute=True) |
---|
| 1628 | #Note the order is unknown, due to using sets |
---|
| 1629 | # and it changes from windows to linux |
---|
| 1630 | try: |
---|
[6153] | 1631 | assert num.allclose(points[1][0], 308728.009) |
---|
| 1632 | assert num.allclose(points[1][1], 6180432.601) |
---|
[6304] | 1633 | assert num.allclose(points[0][0], 222908.705) |
---|
[6153] | 1634 | assert num.allclose(points[0][1], 6233785.284) |
---|
[6086] | 1635 | except AssertionError: |
---|
[6153] | 1636 | assert num.allclose(points[0][0], 308728.009) |
---|
| 1637 | assert num.allclose(points[0][1], 6180432.601) |
---|
[6304] | 1638 | assert num.allclose(points[1][0], 222908.705) |
---|
[6153] | 1639 | assert num.allclose(points[1][1], 6233785.284) |
---|
[6304] | 1640 | |
---|
[6086] | 1641 | self.failUnless(gsd.get_geo_reference().get_zone() == 56, |
---|
| 1642 | 'Bad zone error!') |
---|
| 1643 | points = gsd.get_data_points(as_lat_long=True) |
---|
| 1644 | try: |
---|
[6153] | 1645 | assert num.allclose(points[0][0], -34) |
---|
| 1646 | assert num.allclose(points[0][1], 150) |
---|
[6086] | 1647 | except AssertionError: |
---|
[6153] | 1648 | assert num.allclose(points[1][0], -34) |
---|
| 1649 | assert num.allclose(points[1][1], 150) |
---|
[6086] | 1650 | |
---|
| 1651 | def test_len(self): |
---|
| 1652 | points = [[1.0, 2.1], [3.0, 5.3]] |
---|
| 1653 | G = Geospatial_data(points) |
---|
[6304] | 1654 | self.failUnless(2 == len(G), 'Len error!') |
---|
| 1655 | |
---|
[6086] | 1656 | points = [[1.0, 2.1]] |
---|
| 1657 | G = Geospatial_data(points) |
---|
[6304] | 1658 | self.failUnless(1 == len(G), 'Len error!') |
---|
[6086] | 1659 | |
---|
| 1660 | points = [[1.0, 2.1], [3.0, 5.3], [3.0, 5.3], [3.0, 5.3]] |
---|
| 1661 | G = Geospatial_data(points) |
---|
[6304] | 1662 | self.failUnless(4 == len(G), 'Len error!') |
---|
| 1663 | |
---|
[6086] | 1664 | def test_split(self): |
---|
| 1665 | """test if the results from spilt are disjoin sets""" |
---|
[6304] | 1666 | |
---|
| 1667 | # below is a workaround until randint works on cyclones compute nodes |
---|
| 1668 | if get_host_name()[8:9] != '0': |
---|
[6086] | 1669 | points = [[1.0, 1.0], [1.0, 2.0],[1.0, 3.0], [1.0, 4.0], [1.0, 5.0], |
---|
| 1670 | [2.0, 1.0], [2.0, 2.0],[2.0, 3.0], [2.0, 4.0], [2.0, 5.0], |
---|
| 1671 | [3.0, 1.0], [3.0, 2.0],[3.0, 3.0], [3.0, 4.0], [3.0, 5.0], |
---|
| 1672 | [4.0, 1.0], [4.0, 2.0],[4.0, 3.0], [4.0, 4.0], [4.0, 5.0], |
---|
| 1673 | [5.0, 1.0], [5.0, 2.0],[5.0, 3.0], [5.0, 4.0], [5.0, 5.0]] |
---|
[6304] | 1674 | attributes = {'depth': [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, |
---|
| 1675 | 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, |
---|
| 1676 | 21, 22, 23, 24, 25], |
---|
| 1677 | 'speed': [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, |
---|
| 1678 | 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, |
---|
| 1679 | 21, 22, 23, 24, 25]} |
---|
[6086] | 1680 | G = Geospatial_data(points, attributes) |
---|
[6304] | 1681 | |
---|
[6086] | 1682 | factor = 0.21 |
---|
[6304] | 1683 | |
---|
| 1684 | # will return G1 with 10% of points and G2 with 90% |
---|
| 1685 | G1, G2 = G.split(factor, 100) |
---|
[6153] | 1686 | assert num.allclose(len(G), len(G1)+len(G2)) |
---|
| 1687 | assert num.allclose(round(len(G)*factor), len(G1)) |
---|
[6304] | 1688 | |
---|
[6086] | 1689 | P = G1.get_data_points(absolute=False) |
---|
[6768] | 1690 | expected = [[5.,4.], [4.,1.], [2.,4.], [2.,3.], [1.,4.]] |
---|
| 1691 | # expected = [[5.0, 4.0], [4.0, 3.0], [4.0, 2.0], |
---|
| 1692 | # [3.0, 1.0], [2.0, 3.0]] |
---|
[6304] | 1693 | msg = 'Expected %s, but\nP=%s' % (str(expected), str(P)) |
---|
| 1694 | assert num.allclose(P, expected), msg |
---|
| 1695 | |
---|
[6086] | 1696 | A = G1.get_attributes() |
---|
[6768] | 1697 | expected = [24, 16, 9, 8, 4] |
---|
[6304] | 1698 | msg = 'expected=%s, but A=%s' % (str(expected), str(A)) |
---|
| 1699 | assert num.allclose(A, expected), msg |
---|
| 1700 | |
---|
[6086] | 1701 | def test_split1(self): |
---|
[6304] | 1702 | """test if the results from split are disjoint sets""" |
---|
[6086] | 1703 | |
---|
[6304] | 1704 | # below is a workaround until randint works on cyclones compute nodes |
---|
| 1705 | if get_host_name()[8:9] != '0': |
---|
[6768] | 1706 | from numpy.random import randint, seed |
---|
[6304] | 1707 | |
---|
[6768] | 1708 | seed((100, 100)) |
---|
[6304] | 1709 | a_points = randint(0, 999999, (10,2)) |
---|
[6086] | 1710 | points = a_points.tolist() |
---|
[6304] | 1711 | |
---|
[6086] | 1712 | G = Geospatial_data(points) |
---|
[6304] | 1713 | |
---|
[6086] | 1714 | factor = 0.1 |
---|
[6304] | 1715 | |
---|
| 1716 | # will return G1 with 10% of points and G2 with 90% |
---|
| 1717 | G1, G2 = G.split(factor, 100) |
---|
[6153] | 1718 | assert num.allclose(len(G), len(G1)+len(G2)) |
---|
| 1719 | assert num.allclose(round(len(G)*factor), len(G1)) |
---|
[6304] | 1720 | |
---|
[6086] | 1721 | P = G1.get_data_points(absolute=False) |
---|
[6768] | 1722 | expected = [[425835., 137518.]] |
---|
[6304] | 1723 | msg = 'expected=%s, but\nP=%s' % (str(expected), str(P)) |
---|
| 1724 | assert num.allclose(P, expected), msg |
---|
[6086] | 1725 | |
---|
| 1726 | def test_find_optimal_smoothing_parameter(self): |
---|
| 1727 | """ |
---|
[6304] | 1728 | Creates a elevation file representing a hill (sort of) and runs |
---|
[6086] | 1729 | find_optimal_smoothing_parameter for 3 different alphas, |
---|
[6304] | 1730 | |
---|
[6086] | 1731 | NOTE the random number seed is provided to control the results |
---|
| 1732 | """ |
---|
[6304] | 1733 | |
---|
[6086] | 1734 | from cmath import cos |
---|
| 1735 | |
---|
[6304] | 1736 | # below is a workaround until randint works on cyclones compute nodes |
---|
[6086] | 1737 | if get_host_name()[8:9]!='0': |
---|
[6304] | 1738 | filename = tempfile.mktemp('.csv') |
---|
| 1739 | file = open(filename, 'w') |
---|
| 1740 | file.write('x,y,elevation \n') |
---|
[6086] | 1741 | |
---|
[6304] | 1742 | for i in range(-5, 6): |
---|
| 1743 | for j in range(-5, 6): |
---|
| 1744 | # this equation makes a surface like a circle ripple |
---|
[6086] | 1745 | z = abs(cos(((i*i) + (j*j))*.1)*2) |
---|
[6304] | 1746 | file.write("%s, %s, %s\n" % (i, j, z)) |
---|
| 1747 | |
---|
[6086] | 1748 | file.close() |
---|
[6304] | 1749 | |
---|
| 1750 | value, alpha = find_optimal_smoothing_parameter(data_file=filename, |
---|
[6086] | 1751 | alpha_list=[0.0001, 0.01, 1], |
---|
| 1752 | mesh_file=None, |
---|
| 1753 | mesh_resolution=3, |
---|
| 1754 | north_boundary=5, |
---|
| 1755 | south_boundary=-5, |
---|
| 1756 | east_boundary=5, |
---|
| 1757 | west_boundary=-5, |
---|
| 1758 | plot_name=None, |
---|
| 1759 | seed_num=100000, |
---|
| 1760 | verbose=False) |
---|
[6304] | 1761 | |
---|
[6086] | 1762 | os.remove(filename) |
---|
[6304] | 1763 | |
---|
[6086] | 1764 | # print value, alpha |
---|
[6304] | 1765 | assert(alpha == 0.01) |
---|
[6086] | 1766 | |
---|
| 1767 | def test_find_optimal_smoothing_parameter1(self): |
---|
| 1768 | """ |
---|
[6304] | 1769 | Creates a elevation file representing a hill (sort of) and |
---|
[6086] | 1770 | Then creates a mesh file and passes the mesh file and the elevation |
---|
| 1771 | file to find_optimal_smoothing_parameter for 3 different alphas, |
---|
[6304] | 1772 | |
---|
[6086] | 1773 | NOTE the random number seed is provided to control the results |
---|
| 1774 | """ |
---|
[6304] | 1775 | |
---|
| 1776 | # below is a workaround until randint works on cyclones compute nodes |
---|
[6086] | 1777 | if get_host_name()[8:9]!='0': |
---|
| 1778 | from cmath import cos |
---|
| 1779 | from anuga.pmesh.mesh_interface import create_mesh_from_regions |
---|
[6304] | 1780 | |
---|
| 1781 | filename = tempfile.mktemp('.csv') |
---|
| 1782 | file = open(filename, 'w') |
---|
| 1783 | file.write('x,y,elevation \n') |
---|
| 1784 | |
---|
| 1785 | for i in range(-5 ,6): |
---|
| 1786 | for j in range(-5, 6): |
---|
| 1787 | # this equation makes a surface like a circle ripple |
---|
[6086] | 1788 | z = abs(cos(((i*i) + (j*j))*.1)*2) |
---|
[6304] | 1789 | file.write('%s, %s, %s\n' % (i, j, z)) |
---|
| 1790 | |
---|
[6086] | 1791 | file.close() |
---|
[6304] | 1792 | |
---|
| 1793 | poly=[[5,5], [5,-5], [-5,-5], [-5,5]] |
---|
| 1794 | internal_poly=[[[[1,1], [1,-1], [-1,-1], [-1,1]], .5]] |
---|
| 1795 | mesh_filename= tempfile.mktemp('.msh') |
---|
| 1796 | |
---|
[6086] | 1797 | create_mesh_from_regions(poly, |
---|
[6304] | 1798 | boundary_tags={'back': [2], |
---|
| 1799 | 'side': [1,3], |
---|
| 1800 | 'ocean': [0]}, |
---|
| 1801 | maximum_triangle_area=3, |
---|
| 1802 | interior_regions=internal_poly, |
---|
| 1803 | filename=mesh_filename, |
---|
| 1804 | use_cache=False, |
---|
| 1805 | verbose=False) |
---|
| 1806 | |
---|
| 1807 | value, alpha = find_optimal_smoothing_parameter(data_file=filename, |
---|
[6086] | 1808 | alpha_list=[0.0001, 0.01, 1], |
---|
| 1809 | mesh_file=mesh_filename, |
---|
| 1810 | plot_name=None, |
---|
| 1811 | seed_num=174, |
---|
| 1812 | verbose=False) |
---|
[6304] | 1813 | |
---|
[6086] | 1814 | os.remove(filename) |
---|
| 1815 | os.remove(mesh_filename) |
---|
| 1816 | |
---|
[6768] | 1817 | msg = 'alpha=%s' % str(alpha) |
---|
| 1818 | # 0.01 was expected with Numeric.RandomArray RNG |
---|
| 1819 | assert alpha==1.0, msg |
---|
[6304] | 1820 | |
---|
[6086] | 1821 | def test_find_optimal_smoothing_parameter2(self): |
---|
[6304] | 1822 | '''Tests requirement that mesh file must exist or IOError is thrown |
---|
| 1823 | |
---|
[6086] | 1824 | NOTE the random number seed is provided to control the results |
---|
[6304] | 1825 | ''' |
---|
| 1826 | |
---|
[6086] | 1827 | from cmath import cos |
---|
| 1828 | from anuga.pmesh.mesh_interface import create_mesh_from_regions |
---|
[6304] | 1829 | |
---|
| 1830 | filename = tempfile.mktemp('.csv') |
---|
| 1831 | mesh_filename= tempfile.mktemp('.msh') |
---|
| 1832 | |
---|
[6086] | 1833 | try: |
---|
[6304] | 1834 | value, alpha = find_optimal_smoothing_parameter(data_file=filename, |
---|
[6086] | 1835 | alpha_list=[0.0001, 0.01, 1], |
---|
| 1836 | mesh_file=mesh_filename, |
---|
| 1837 | plot_name=None, |
---|
| 1838 | seed_num=174, |
---|
| 1839 | verbose=False) |
---|
| 1840 | except IOError: |
---|
| 1841 | pass |
---|
| 1842 | else: |
---|
[6304] | 1843 | self.fail('Error not thrown error!') |
---|
| 1844 | |
---|
| 1845 | ################################################################################ |
---|
| 1846 | |
---|
[6086] | 1847 | if __name__ == "__main__": |
---|
| 1848 | suite = unittest.makeSuite(Test_Geospatial_data, 'test') |
---|
| 1849 | runner = unittest.TextTestRunner() #verbosity=2) |
---|
| 1850 | runner.run(suite) |
---|
| 1851 | |
---|
[6304] | 1852 | |
---|