Context Navigation

source: anuga_core/source_numpy_conversion/anuga/geospatial_data/test_geospatial_data.py @ 5915

Last change on this file since 5915 was 5915, checked in by rwilson, 15 years ago
NumPy? conversion.
File size: 73.2 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: