Context Navigation

source: branches/source_numpy_conversion/anuga/geospatial_data/test_geospatial_data.py @ 6768

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

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

Download in other formats: