Context Navigation

source: anuga_core/source/anuga/geospatial_data/geospatial_data.py @ 4095

Last change on this file since 4095 was 4095, checked in by nick, 17 years ago
add number of trigs in mesh
File size: 35.6 KB

Line
1	"""Class Geospatial_data - Manipulation of locations on the planet and
2	associated attributes.
3
4	"""
5
6	from os import access, F_OK, R_OK
7	from types import DictType
8
9	from Numeric import concatenate, array, Float, shape, reshape, ravel, take, \
10	size, shape
11	from random import randint
12	#from MA import tolist
13
14	from anuga.utilities.numerical_tools import ensure_numeric
15	from anuga.coordinate_transforms.geo_reference import Geo_reference, TitleError
16	from anuga.coordinate_transforms.redfearn import convert_from_latlon_to_utm
17
18
19	class Geospatial_data:
20
21	def __init__(self,
22	data_points = None,
23	attributes = None,
24	geo_reference = None,
25	default_attribute_name = None,
26	file_name = None,
27	delimiter = None,
28	latitudes = None,
29	longitudes = None,
30	points_are_lats_longs = False,
31	verbose = False):
32
33
34	"""
35	Create instance from data points and associated attributes
36
37	data_points: x,y coordinates in meters. Type must be either a
38	sequence of 2-tuples or an Mx2 Numeric array of floats.
39
40	attributes: Associated values for each data point. The type
41	must be either a list or an array of length M or a dictionary
42	of lists (or arrays) of length M. In the latter case the keys
43	in the dictionary represent the attribute names, in the former
44	the attribute will get the default name "attribute".
45
46	geo_reference: Object representing the origin of the data
47	points. It contains UTM zone, easting and northing and data
48	points are assumed to be relative to this origin.
49	If geo_reference is None, the default geo ref object is used
50
51	default_attribute_name: Name of default attribute to be used with
52	get_attribute_values. The idea is that the dataset can be
53	equipped with information about which attribute to return.
54	If None, the default is the "first"
55
56	file_name: Name of input netCDF file or xya file. netCDF file must
57	have dimensions "points" etc.
58	xya file is a comma seperated file with x, y and attribute data.
59	the first line must be the attribute names eg elevation
60
61	The format for a .xya file is:
62	1st line: [attribute names]
63	other lines: x y [attributes]
64
65	for example:
66	elevation, friction
67	0.6, 0.7, 4.9, 0.3
68	1.9, 2.8, 5, 0.3
69	2.7, 2.4, 5.2, 0.3
70
71	The first two columns are always implicitly assumed to be x, y coordinates.
72	Use the same delimiter for the attribute names and the data
73
74	An xya file can optionally end with
75	#geo reference
76	56
77	466600.0
78	8644444.0
79
80	When the 1st # is the zone,
81	2nd # the xllcorner and
82	3rd # the yllcorner
83
84	An issue with the xya format is that the attribute column order
85	is not be controlled. The info is stored in a dictionary and it's
86	written however
87
88	The format for a Points dictionary is:
89
90	['pointlist'] a 2 column array describing points. 1st column x, 2nd column y.
91	['attributelist'], a dictionary of 1D arrays, representing attribute values
92	at the point. The dictionary key is the attribute header.
93	['geo_reference'] a Geo_refernece object. Use if the point information
94	is relative. This is optional.
95	eg
96	dic['pointlist'] = [[1.0,2.0],[3.0,5.0]]
97	dic['attributelist']['elevation'] = [[7.0,5.0]
98
99	delimiter: is the file delimiter that will be used when
100	importing the file
101
102	verbose:
103
104	"""
105
106	if isinstance(data_points, basestring):
107	# assume data point is really a file name
108	file_name = data_points
109
110	self.set_verbose(verbose)
111	self.geo_reference=None #create the attribute
112	if file_name is None:
113	if delimiter is not None:
114	msg = 'No file specified yet a delimiter is provided!'
115	raise ValueError, msg
116	file_name = None
117	if latitudes is not None or longitudes is not None or \
118	points_are_lats_longs:
119	data_points, geo_reference = \
120	self._set_using_lat_long(latitudes=latitudes,
121	longitudes=longitudes,
122	geo_reference=geo_reference,
123	data_points=data_points,
124	points_are_lats_longs=points_are_lats_longs)
125	self.check_data_points(data_points)
126	self.set_attributes(attributes)
127	self.set_geo_reference(geo_reference)
128	self.set_default_attribute_name(default_attribute_name)
129
130	else:
131	# watch for case where file name and points, attributes etc are provided!!
132	# if file name then all provided info will be removed!
133	self.import_points_file(file_name, delimiter, verbose)
134
135	self.check_data_points(self.data_points)
136	self.set_attributes(self.attributes)
137	self.set_geo_reference(self.geo_reference)
138	self.set_default_attribute_name(default_attribute_name)
139
140
141	assert self.attributes is None or isinstance(self.attributes, DictType)
142
143
144	def __len__(self):
145	return len(self.data_points)
146
147	def __repr__(self):
148	return str(self.get_data_points(absolute=True))
149
150
151	def check_data_points(self, data_points):
152	"""Checks data points
153	"""
154
155	if data_points is None:
156	self.data_points = None
157	msg = 'There is no data or file provided!'
158	raise ValueError, msg
159
160	else:
161	self.data_points = ensure_numeric(data_points)
162	assert len(self.data_points.shape) == 2
163	assert self.data_points.shape[1] == 2
164
165	def set_attributes(self, attributes):
166	"""Check and assign attributes dictionary
167	"""
168
169	if attributes is None:
170	self.attributes = None
171	return
172
173	if not isinstance(attributes, DictType):
174	#Convert single attribute into dictionary
175	attributes = {'attribute': attributes}
176
177	#Check input attributes
178	for key in attributes.keys():
179	try:
180	attributes[key] = ensure_numeric(attributes[key])
181	except:
182	msg = 'Attribute %s could not be converted' %key
183	msg += 'to a numeric vector'
184	raise msg
185
186	self.attributes = attributes
187
188
189	def set_geo_reference(self, geo_reference):
190
191	from anuga.coordinate_transforms.geo_reference import Geo_reference
192
193	if geo_reference is None:
194	geo_reference = Geo_reference() # Use default
195	if not isinstance(geo_reference, Geo_reference):
196	msg = 'Argument geo_reference must be a valid Geo_reference \n'
197	msg += 'object or None.'
198	raise msg
199
200	# if a geo ref already exists, change the point data to
201	# represent the new geo-ref
202	if self.geo_reference is not None:
203	#FIXME: Maybe put out a warning here...
204	self.data_points = self.get_data_points \
205	(geo_reference=geo_reference)
206
207	self.geo_reference = geo_reference
208
209
210	def set_default_attribute_name(self, default_attribute_name):
211	self.default_attribute_name = default_attribute_name
212
213	def set_verbose(self, verbose = False):
214	if verbose is not False:
215	verbose = True
216	else:
217	verbose = False
218
219	def clip(self, polygon, closed=True):
220	"""Clip geospatial data by a polygon
221
222	Input
223	polygon - Either a list of points, an Nx2 array or
224	a Geospatial data object.
225	closed - (optional) determine whether points on boundary should be
226	regarded as belonging to the polygon (closed = True)
227	or not (closed = False). Default is True.
228
229	Output
230	New geospatial data object representing points inside
231	specified polygon.
232
233	"""
234
235	from anuga.utilities.polygon import inside_polygon
236
237	if isinstance(polygon, Geospatial_data):
238	# Polygon is an object - extract points
239	polygon = polygon.get_data_points()
240
241	points = self.get_data_points()
242	inside_indices = inside_polygon(points, polygon, closed)
243
244	clipped_G = self.get_sample(inside_indices)
245	# clipped_points = take(points, inside_indices)
246
247	# Clip all attributes
248	# attributes = self.get_all_attributes()
249
250	# clipped_attributes = {}
251	# if attributes is not None:
252	# for key, att in attributes.items():
253	# clipped_attributes[key] = take(att, inside_indices)
254
255	# return Geospatial_data(clipped_points, clipped_attributes)
256	return clipped_G
257
258
259	def clip_outside(self, polygon, closed=True):
260	"""Clip geospatial date by a polygon, keeping data OUTSIDE of polygon
261
262	Input
263	polygon - Either a list of points, an Nx2 array or
264	a Geospatial data object.
265	closed - (optional) determine whether points on boundary should be
266	regarded as belonging to the polygon (closed = True)
267	or not (closed = False). Default is True.
268
269	Output
270	Geospatial data object representing point OUTSIDE specified polygon
271
272	"""
273
274	from anuga.utilities.polygon import outside_polygon
275
276	if isinstance(polygon, Geospatial_data):
277	# Polygon is an object - extract points
278	polygon = polygon.get_data_points()
279
280	points = self.get_data_points()
281	outside_indices = outside_polygon(points, polygon, closed)
282
283	clipped_G = self.get_sample(outside_indices)
284
285	# clipped_points = take(points, outside_indices)
286
287	# Clip all attributes
288	# attributes = self.get_all_attributes()
289
290	# clipped_attributes = {}
291	# if attributes is not None:
292	# for key, att in attributes.items():
293	# clipped_attributes[key] = take(att, outside_indices)
294
295	# return Geospatial_data(clipped_points, clipped_attributes)
296	return clipped_G
297
298
299	def _set_using_lat_long(self,
300	latitudes,
301	longitudes,
302	geo_reference,
303	data_points,
304	points_are_lats_longs):
305
306	if geo_reference is not None:
307	msg = """A georeference is specified yet latitude and longitude
308	are also specified!"""
309	raise ValueError, msg
310
311	if data_points is not None and not points_are_lats_longs:
312	msg = """Data points are specified yet latitude and
313	longitude are also specified!"""
314	raise ValueError, msg
315
316	if points_are_lats_longs:
317	if data_points is None:
318	msg = """Data points are not specified !"""
319	raise ValueError, msg
320	lats_longs = ensure_numeric(data_points)
321	latitudes = ravel(lats_longs[:,0:1])
322	longitudes = ravel(lats_longs[:,1:])
323
324	if latitudes is None and longitudes is None:
325	msg = """Latitudes and Longitudes are not."""
326	raise ValueError, msg
327
328	if latitudes is None:
329	msg = """Longitudes are specified yet latitudes aren't!"""
330	raise ValueError, msg
331
332	if longitudes is None:
333	msg = """Latitudes are specified yet longitudes aren't!"""
334	raise ValueError, msg
335
336	data_points, zone = convert_from_latlon_to_utm(latitudes=latitudes,
337	longitudes=longitudes)
338
339	return data_points, Geo_reference(zone=zone)
340
341	def get_geo_reference(self):
342	return self.geo_reference
343
344	def get_data_points(self, absolute = True, geo_reference=None):
345	"""Get coordinates for all data points as an Nx2 array
346
347	If absolute is True absolute UTM coordinates are returned otherwise
348	returned coordinates are relative to the internal georeference's
349	xll and yll corners.
350
351	If a geo_reference is passed the points are returned relative
352	to that geo_reference.
353
354	Default: absolute is True.
355	"""
356
357	if absolute is True and geo_reference is None:
358	return self.geo_reference.get_absolute(self.data_points)
359	elif geo_reference is not None:
360	return geo_reference.change_points_geo_ref \
361	(self.data_points,
362	self.geo_reference)
363	else:
364	return self.data_points
365
366
367	def get_attributes(self, attribute_name = None):
368	"""Return values for one named attribute.
369
370	If attribute_name is None, default_attribute_name is used
371	"""
372
373	if attribute_name is None:
374	if self.default_attribute_name is not None:
375	attribute_name = self.default_attribute_name
376	else:
377	attribute_name = self.attributes.keys()[0]
378	# above line takes the first one from keys
379
380
381	msg = 'Attribute name %s does not exist in data set' %attribute_name
382	assert self.attributes.has_key(attribute_name), msg
383
384	return self.attributes[attribute_name]
385
386	def get_all_attributes(self):
387	"""
388	Return values for all attributes.
389	The return value is either None or a dictionary (possibly empty).
390	"""
391
392	return self.attributes
393
394	def __add__(self, other):
395	"""
396	Returns the addition of 2 geospatical objects,
397	objects are concatencated to the end of each other
398
399	NOTE: doesn't add if objects contain different
400	attributes
401
402	Always return relative points!
403	"""
404
405	# find objects zone and checks if the same
406	geo_ref1 = self.get_geo_reference()
407	zone1 = geo_ref1.get_zone()
408
409	geo_ref2 = other.get_geo_reference()
410	zone2 = geo_ref2.get_zone()
411
412	geo_ref1.reconcile_zones(geo_ref2)
413
414
415	# sets xll and yll as the smallest from self and other
416	# FIXME (Duncan and Ole): use lower left corner derived from
417	# absolute coordinates
418	if self.geo_reference.xllcorner <= other.geo_reference.xllcorner:
419	xll = self.geo_reference.xllcorner
420	else:
421	xll = other.geo_reference.xllcorner
422
423	if self.geo_reference.yllcorner <= other.geo_reference.yllcorner:
424	yll = self.geo_reference.yllcorner
425	else:
426	yll = other.geo_reference.yllcorner
427	new_geo_ref = Geo_reference(geo_ref1.get_zone(), xll, yll)
428
429	xll = yll = 0.
430
431	relative_points1 = self.get_data_points(absolute = False)
432	relative_points2 = other.get_data_points(absolute = False)
433
434
435	new_relative_points1 = new_geo_ref.\
436	change_points_geo_ref(relative_points1,
437	geo_ref1)
438	new_relative_points2 = new_geo_ref.\
439	change_points_geo_ref(relative_points2,
440	geo_ref2)
441
442	# Now both point sets are relative to new_geo_ref and
443	# zones have been reconciled
444
445	# Concatenate points
446	new_points = concatenate((new_relative_points1,
447	new_relative_points2),
448	axis = 0)
449
450	# Concatenate attributes if any
451	if self.attributes is None:
452	if other.attributes is not None:
453	msg = 'Both geospatial_data objects must have the same \n'
454	msg += 'attributes to allow addition.'
455	raise Exception, msg
456
457	new_attributes = None
458	else:
459	new_attributes = {}
460	for x in self.attributes.keys():
461	if other.attributes.has_key(x):
462
463	attrib1 = self.attributes[x]
464	attrib2 = other.attributes[x]
465	new_attributes[x] = concatenate((attrib1, attrib2))
466
467	else:
468	msg = 'Both geospatial_data objects must have the same \n'
469	msg += 'attributes to allow addition.'
470	raise Exception, msg
471
472	# Instantiate new data object and return
473	return Geospatial_data(new_points,
474	new_attributes,
475	new_geo_ref)
476
477	###
478	# IMPORT/EXPORT POINTS FILES
479	###
480
481	def import_points_file(self, file_name, delimiter = None, verbose = False):
482	""" load an .xya or .pts file
483	Note: will throw an IOError if it can't load the file.
484	Catch these!
485
486	Post condition: self.attributes dictionary has been set
487	"""
488
489	if access(file_name, F_OK) == 0 :
490	msg = 'File %s does not exist or is not accessible' %file_name
491	raise IOError, msg
492
493	attributes = {}
494	if file_name[-4:]== ".xya":
495	try:
496	if delimiter == None:
497	try:
498	fd = open(file_name)
499	data_points, attributes, geo_reference = _read_xya_file(fd, ',')
500	except TitleError:
501	fd.close()
502	fd = open(file_name)
503	data_points, attributes, geo_reference = _read_xya_file(fd, ' ')
504	else:
505	fd = open(file_name)
506	data_points, attributes, geo_reference = _read_xya_file(fd, delimiter)
507	fd.close()
508	except (IndexError,ValueError,SyntaxError):
509	fd.close()
510	msg = 'Could not open file %s ' %file_name
511	raise IOError, msg
512	except IOError, e:
513	fd.close()
514	# Catch this to add an error message
515	msg = 'Could not open file or incorrect file format %s:%s' %(file_name, e)
516	raise IOError, msg
517
518	elif file_name[-4:]== ".pts":
519	try:
520	data_points, attributes, geo_reference = _read_pts_file(file_name, verbose)
521	except IOError, e:
522	msg = 'Could not open file %s ' %file_name
523	raise IOError, msg
524
525	elif file_name[-4:]== ".xxx":
526	#let's do ticket#116 stuff
527	#
528	try:
529	data_points, attributes, geo_reference = _read_csv_file(file_name, verbose)
530	except IOError, e:
531	msg = 'Could not open file %s ' %file_name
532	raise IOError, msg
533	else:
534	msg = 'Extension %s is unknown' %file_name[-4:]
535	raise IOError, msg
536	# print'in import data_points', data_points
537	# print'in import attributes', attributes
538	# print'in import data_points', geo_reference
539	self.data_points = data_points
540	self.attributes = attributes
541	self.geo_reference = geo_reference
542
543	# return all_data
544
545	def export_points_file(self, file_name, absolute=True):
546
547	"""
548	write a points file, file_name, as a text (.xya) or binary (.pts) file
549	file_name is the file name, including the extension
550	The point_dict is defined at the top of this file.
551
552	If absolute is True data points at returned added to the xll and yll
553	and geo_reference as None
554
555	If absolute is False data points at returned as relative to the xll
556	and yll and geo_reference remains uneffected
557	"""
558
559	if (file_name[-4:] == ".xya"):
560	if absolute is True:
561	_write_xya_file(file_name,
562	self.get_data_points(absolute=True),
563	self.get_all_attributes())
564	else:
565	_write_xya_file(file_name,
566	self.get_data_points(absolute=False),
567	self.get_all_attributes(),
568	self.get_geo_reference())
569
570	elif (file_name[-4:] == ".pts"):
571	if absolute is True:
572	_write_pts_file(file_name,
573	self.get_data_points(absolute),
574	self.get_all_attributes())
575	else:
576	_write_pts_file(file_name,
577	self.get_data_points(absolute),
578	self.get_all_attributes(),
579	self.get_geo_reference())
580	else:
581	msg = 'Unknown file type %s ' %file_name
582	raise IOError, msg
583
584	def get_sample(self, indices):
585	""" Returns a object which is a subset of the original
586	and the data points and attributes in this new object refer to
587	the indices provided
588
589	Input
590	indices- a list of integers that represent the new object
591	Output
592	New geospatial data object representing points specified by
593	the indices
594	"""
595	#FIXME: add the geo_reference to this
596
597	points = self.get_data_points()
598	sampled_points = take(points, indices)
599
600	attributes = self.get_all_attributes()
601
602	sampled_attributes = {}
603	if attributes is not None:
604	for key, att in attributes.items():
605	sampled_attributes[key] = take(att, indices)
606
607	return Geospatial_data(sampled_points, sampled_attributes)
608
609
610	def split(self, factor=0.5):
611	"""Returns two geospatial_data object, first is size of the 'factor'
612	smaller the original and the second is the remainer. The two new
613	object are disjoin set of each other.
614
615	Points of the two new object have selected RANDOMLY.
616	AND if factor is a decimal it will round (2.25 to 2 and 2.5 to 3)
617
618
619	Input - the factor which to split the object, if 0.1 then 10% of the
620	object will be returned
621
622	Output - two geospatial_data objects that are disjoint sets of the
623	original
624	"""
625
626	i=0
627	self_size = len(self)
628	random_list = []
629	remainder_list = []
630	new_size = round(factor*self_size)
631	# print'Split original %s by %s' %(self_size, factor)
632	# print'New samples are %s and %s in size' %(int(round(factor*self_size)),int(self_size-new_size))
633
634	#find unique random numbers
635	while i < new_size:
636	random_num = randint(0,self_size-1)
637	if random_num not in random_list:
638	random_list.append(random_num)
639	i=i+1
640
641	#Make list of opposite to random_list
642	for i in range(0,self_size,1):
643	remainder_list.append(i)
644
645	#remove random list from remainder_list to get correct remainder_list
646	#need to sort and reverse so the pop() works correctly
647	random_list.sort()
648	random_list.reverse()
649	for i in random_list:
650	remainder_list.pop(i)
651
652	#get new samples
653	G1 = self.get_sample(random_list)
654	G2 = self.get_sample(remainder_list)
655
656	return G1, G2
657
658
659	def _read_pts_file(file_name, verbose = False):
660	"""Read .pts NetCDF file
661
662	Return a dic of array of points, and dic of array of attribute
663	eg
664	dic['points'] = [[1.0,2.0],[3.0,5.0]]
665	dic['attributelist']['elevation'] = [[7.0,5.0]
666	"""
667
668	from Scientific.IO.NetCDF import NetCDFFile
669
670	if verbose: print 'Reading ', file_name
671
672
673	# see if the file is there. Throw a QUIET IO error if it isn't
674	fd = open(file_name,'r')
675	fd.close()
676
677	#throws prints to screen if file not present
678	fid = NetCDFFile(file_name, 'r')
679
680	# point_atts = {}
681	# Get the variables
682	# point_atts['pointlist'] = array(fid.variables['points'])
683	pointlist = array(fid.variables['points'])
684	keys = fid.variables.keys()
685	if verbose: print 'Got %d variables: %s' %(len(keys), keys)
686	try:
687	keys.remove('points')
688	except IOError, e:
689	fid.close()
690	msg = 'Expected keyword "points" but could not find it'
691	raise IOError, msg
692
693	attributes = {}
694	for key in keys:
695	if verbose: print "reading attribute '%s'" %key
696
697	attributes[key] = array(fid.variables[key])
698
699	# point_atts['attributelist'] = attributes
700
701	try:
702	geo_reference = Geo_reference(NetCDFObject=fid)
703	# point_atts['geo_reference'] = geo_reference
704	except AttributeError, e:
705	#geo_ref not compulsory
706	# point_atts['geo_reference'] = None
707	geo_reference = None
708
709	fid.close()
710
711	return pointlist, attributes, geo_reference
712
713
714	def _read_csv_file(file_name, verbose = False):
715	"""Read .csv file
716
717	Return a dic of array of points, and dic of array of attribute
718	eg
719	dic['points'] = [[1.0,2.0],[3.0,5.0]]
720	dic['attributelist']['elevation'] = [[7.0,5.0]
721	"""
722
723	from anuga.shallow_water.data_manager import Exposure_csv
724	csv =Exposure_csv(file_name)
725
726	return pointlist, attributes, geo_reference
727
728	def _read_xya_file( fd, delimiter):
729	points = []
730	pointattributes = []
731	title = fd.readline()
732	att_names = clean_line(title,delimiter)
733	att_dict = {}
734	line = fd.readline()
735	numbers = clean_line(line,delimiter)
736
737	while len(numbers) > 1 and line[0] <> '#':
738	if numbers != []:
739	try:
740	x = float(numbers[0])
741	y = float(numbers[1])
742	points.append([x,y])
743	numbers.pop(0)
744	numbers.pop(0)
745	if len(att_names) != len(numbers):
746	fd.close()
747	# It might not be a problem with the title
748	#raise TitleAmountError
749	raise IOError
750	for i,num in enumerate(numbers):
751	num.strip()
752	if num != '\n' and num != '':
753	#attributes.append(float(num))
754	att_dict.setdefault(att_names[i],[]).append(float(num))
755	except ValueError:
756	raise SyntaxError
757	line = fd.readline()
758	numbers = clean_line(line,delimiter)
759
760	if line == '':
761	geo_reference = None
762	else:
763	geo_reference = Geo_reference(ASCIIFile=fd,read_title=line)
764
765
766	pointlist = array(points).astype(Float)
767	for key in att_dict.keys():
768	att_dict[key] = array(att_dict[key]).astype(Float)
769
770	return pointlist, att_dict, geo_reference
771
772	def _write_pts_file(file_name,
773	write_data_points,
774	write_attributes=None,
775	write_geo_reference=None):
776	"""
777	Write .pts NetCDF file
778
779	NOTE: Below might not be valid ask Duncan : NB 5/2006
780
781	WARNING: This function mangles the point_atts data structure
782	#F??ME: (DSG)This format has issues.
783	# There can't be an attribute called points
784	# consider format change
785	# method changed by NB not sure if above statement is correct
786
787	should create new test for this
788	legal_keys = ['pointlist', 'attributelist', 'geo_reference']
789	for key in point_atts.keys():
790	msg = 'Key %s is illegal. Valid keys are %s' %(key, legal_keys)
791	assert key in legal_keys, msg
792	"""
793	from Scientific.IO.NetCDF import NetCDFFile
794	# NetCDF file definition
795	outfile = NetCDFFile(file_name, 'w')
796
797	#Create new file
798	outfile.institution = 'Geoscience Australia'
799	outfile.description = 'NetCDF format for compact and portable storage ' +\
800	'of spatial point data'
801
802	# dimension definitions
803	shape = write_data_points.shape[0]
804	outfile.createDimension('number_of_points', shape)
805	outfile.createDimension('number_of_dimensions', 2) #This is 2d data
806
807	# variable definition
808	outfile.createVariable('points', Float, ('number_of_points',
809	'number_of_dimensions'))
810
811	#create variables
812	outfile.variables['points'][:] = write_data_points #.astype(Float32)
813
814	if write_attributes is not None:
815	for key in write_attributes.keys():
816	outfile.createVariable(key, Float, ('number_of_points',))
817	outfile.variables[key][:] = write_attributes[key] #.astype(Float32)
818
819	if write_geo_reference is not None:
820	write_geo_reference.write_NetCDF(outfile)
821
822	outfile.close()
823
824
825
826	def _write_xya_file(file_name,
827	write_data_points,
828	write_attributes=None,
829	write_geo_reference=None,
830	delimiter = ','):
831	"""
832	export a file, file_name, with the xya format
833
834	"""
835	points = write_data_points
836	pointattributes = write_attributes
837
838	fd = open(file_name,'w')
839	titlelist = ""
840	if pointattributes is not None:
841	for title in pointattributes.keys():
842	titlelist = titlelist + title + delimiter
843	titlelist = titlelist[0:-len(delimiter)] # remove the last delimiter
844	fd.write(titlelist+"\n")
845
846	#<vertex #> <x> <y> [attributes]
847	for i, vert in enumerate( points):
848
849
850	if pointattributes is not None:
851	attlist = ","
852	for att in pointattributes.keys():
853	attlist = attlist + str(pointattributes[att][i])+ delimiter
854	attlist = attlist[0:-len(delimiter)] # remove the last delimiter
855	attlist.strip()
856	else:
857	attlist = ''
858
859	fd.write(str(vert[0]) + delimiter +
860	str(vert[1]) + attlist + "\n")
861
862	if write_geo_reference is not None:
863	write_geo_reference.write_ASCII(fd)
864	fd.close()
865
866
867
868	def _point_atts2array(point_atts):
869	point_atts['pointlist'] = array(point_atts['pointlist']).astype(Float)
870
871	for key in point_atts['attributelist'].keys():
872	point_atts['attributelist'][key]= array(point_atts['attributelist'][key]).astype(Float)
873	return point_atts
874
875
876
877
878	def geospatial_data2points_dictionary(geospatial_data):
879	"""Convert geospatial data to points_dictionary
880	"""
881
882	points_dictionary = {}
883	points_dictionary['pointlist'] = geospatial_data.data_points
884
885	points_dictionary['attributelist'] = {}
886
887	for attribute_name in geospatial_data.attributes.keys():
888	val = geospatial_data.attributes[attribute_name]
889	points_dictionary['attributelist'][attribute_name] = val
890
891	points_dictionary['geo_reference'] = geospatial_data.geo_reference
892
893	return points_dictionary
894
895
896	def points_dictionary2geospatial_data(points_dictionary):
897	"""Convert points_dictionary to geospatial data object
898	"""
899
900	msg = 'Points dictionary must have key pointlist'
901	assert points_dictionary.has_key('pointlist'), msg
902
903	msg = 'Points dictionary must have key attributelist'
904	assert points_dictionary.has_key('attributelist'), msg
905
906	if points_dictionary.has_key('geo_reference'):
907	geo = points_dictionary['geo_reference']
908	else:
909	geo = None
910
911	return Geospatial_data(points_dictionary['pointlist'],
912	points_dictionary['attributelist'],
913	geo_reference = geo)
914
915	def clean_line(line,delimiter):
916	"""Remove whitespace
917	"""
918	#print ">%s" %line
919	line = line.strip()
920	#print "stripped>%s" %line
921	numbers = line.split(delimiter)
922	i = len(numbers) - 1
923	while i >= 0:
924	if numbers[i] == '':
925	numbers.pop(i)
926	else:
927	numbers[i] = numbers[i].strip()
928
929	i += -1
930	#for num in numbers:
931	# print "num>%s<" %num
932	return numbers
933
934	def ensure_absolute(points, geo_reference = None):
935	"""
936	This function inputs several formats and
937	outputs one format. - a numeric array of absolute points.
938
939	Inputed formats are;
940	points: List or numeric array of coordinate pairs [xi, eta] of
941	points or geospatial object or points file name
942
943	mesh_origin: A geo_reference object or 3-tuples consisting of
944	UTM zone, easting and northing.
945	If specified vertex coordinates are assumed to be
946	relative to their respective origins.
947	"""
948	if isinstance(points,type('')):
949	#It's a string
950	#assume it is a point file
951	points = Geospatial_data(file_name = points)
952
953	if isinstance(points,Geospatial_data):
954	points = points.get_data_points( \
955	absolute = True)
956	msg = "Use a Geospatial_data object or a mesh origin. Not both."
957	assert geo_reference == None, msg
958
959	else:
960	points = ensure_numeric(points, Float)
961	if geo_reference is None:
962	geo = None #Geo_reference()
963	else:
964	if isinstance(geo_reference, Geo_reference):
965	geo = geo_reference
966	else:
967	geo = Geo_reference(geo_reference[0],
968	geo_reference[1],
969	geo_reference[2])
970	points = geo.get_absolute(points)
971	return points
972
973
974	def ensure_geospatial(points, geo_reference = None):
975	"""
976	This function inputs several formats and
977	outputs one format. - a geospatial_data instance.
978
979	Inputed formats are;
980	points: List or numeric array of coordinate pairs [xi, eta] of
981	points or geospatial object
982
983	mesh_origin: A geo_reference object or 3-tuples consisting of
984	UTM zone, easting and northing.
985	If specified vertex coordinates are assumed to be
986	relative to their respective origins.
987	"""
988	if isinstance(points,Geospatial_data):
989	msg = "Use a Geospatial_data object or a mesh origin. Not both."
990	assert geo_reference == None, msg
991	return points
992	else:
993	points = ensure_numeric(points, Float)
994	if geo_reference is None:
995	geo = None
996	else:
997	if isinstance(geo_reference, Geo_reference):
998	geo = geo_reference
999	else:
1000	geo = Geo_reference(geo_reference[0],
1001	geo_reference[1],
1002	geo_reference[2])
1003	points = Geospatial_data(data_points=points, geo_reference=geo)
1004	return points
1005
1006	#def file2xya(filename):
1007
1008	# G = Geospatial_data(filename)
1009	# G.export_points_file()
1010
1011

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

Download in other formats: