Context Navigation

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

Last change on this file since 3969 was 3969, checked in by nick, 17 years ago
added split and get sample modules to geospatial_data also removed xxx_add_data_points this was not being used and is old
File size: 34.8 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 RandomArray 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	else:
525	msg = 'Extension %s is unknown' %file_name[-4:]
526	raise IOError, msg
527
528	# print'in import data_points', data_points
529	# print'in import attributes', attributes
530	# print'in import data_points', geo_reference
531	self.data_points = data_points
532	self.attributes = attributes
533	self.geo_reference = geo_reference
534
535	# return all_data
536
537	def export_points_file(self, file_name, absolute=True):
538
539	"""
540	write a points file, file_name, as a text (.xya) or binary (.pts) file
541	file_name is the file name, including the extension
542	The point_dict is defined at the top of this file.
543
544	If absolute is True data points at returned added to the xll and yll
545	and geo_reference as None
546
547	If absolute is False data points at returned as relative to the xll
548	and yll and geo_reference remains uneffected
549	"""
550
551	if (file_name[-4:] == ".xya"):
552	if absolute is True:
553	_write_xya_file(file_name,
554	self.get_data_points(absolute=True),
555	self.get_all_attributes())
556	else:
557	_write_xya_file(file_name,
558	self.get_data_points(absolute=False),
559	self.get_all_attributes(),
560	self.get_geo_reference())
561
562	elif (file_name[-4:] == ".pts"):
563	if absolute is True:
564	_write_pts_file(file_name,
565	self.get_data_points(absolute),
566	self.get_all_attributes())
567	else:
568	_write_pts_file(file_name,
569	self.get_data_points(absolute),
570	self.get_all_attributes(),
571	self.get_geo_reference())
572	else:
573	msg = 'Unknown file type %s ' %file_name
574	raise IOError, msg
575
576	def get_sample(self, indices):
577	""" Returns a object which is a subset of the original
578	and the data points and attributes in this new object refer to
579	the indices provided
580
581	Input
582	indices- a list of integers that represent the new object
583	Output
584	New geospatial data object representing points specified by
585	the indices
586	"""
587	#FIXME: add the geo_reference to this
588
589	points = self.get_data_points()
590	sampled_points = take(points, indices)
591
592	attributes = self.get_all_attributes()
593
594	sampled_attributes = {}
595	if attributes is not None:
596	for key, att in attributes.items():
597	sampled_attributes[key] = take(att, indices)
598
599	return Geospatial_data(sampled_points, sampled_attributes)
600
601
602	def split(self, factor=0.5):
603	"""Returns two geospatial_data object, first is size of the 'factor'
604	smaller the original one and the second is the remainer. The two new
605	object are disjoin set of each other.
606
607	Points of the two new object have selected RANDOMLY.
608	AND if factor is a decimal it will round (2.25 to 2 and 2.5 to 3)
609
610
611	Input - the factor which to split the object, if 0.1 then 10% of the
612	object will be returned
613
614	Output - two geospatial_data objects that are disjoint sets of the
615	original
616	"""
617
618	i=0
619	self_size = len(self)
620	# print 'size', self_size
621	random_list = []
622	remainder_list = []
623	new_size = round(factor*self_size)
624	# print'Split original %s by %s' %(self_size, factor)
625	# print'New samples are %s and %s in size' %(int(round(factor*self_size)),int(self_size-new_size))
626
627	#find unique random numbers
628	while i < new_size:
629	random_num = randint(0,self_size)
630	if random_num not in random_list:
631	random_list.append(random_num)
632	i=i+1
633
634	#Make list of opposite to random_list
635	for i in range(0,self_size,1):
636	remainder_list.append(i)
637	# print 'start remainer',remainder_list
638
639	#need to sort and reverse so the pop() works correctly
640	random_list.sort()
641	random_list.reverse()
642	for i in random_list:
643	remainder_list.pop(i)
644
645	#get new samples
646	G1 = self.get_sample(random_list)
647	G2 = self.get_sample(remainder_list)
648
649	return G1, G2
650
651
652	def _read_pts_file(file_name, verbose = False):
653	"""Read .pts NetCDF file
654
655	Return a dic of array of points, and dic of array of attribute
656	eg
657	dic['points'] = [[1.0,2.0],[3.0,5.0]]
658	dic['attributelist']['elevation'] = [[7.0,5.0]
659	"""
660
661	from Scientific.IO.NetCDF import NetCDFFile
662
663	if verbose: print 'Reading ', file_name
664
665
666	# see if the file is there. Throw a QUIET IO error if it isn't
667	fd = open(file_name,'r')
668	fd.close()
669
670	#throws prints to screen if file not present
671	fid = NetCDFFile(file_name, 'r')
672
673	# point_atts = {}
674	# Get the variables
675	# point_atts['pointlist'] = array(fid.variables['points'])
676	pointlist = array(fid.variables['points'])
677	keys = fid.variables.keys()
678	if verbose: print 'Got %d variables: %s' %(len(keys), keys)
679	try:
680	keys.remove('points')
681	except IOError, e:
682	fid.close()
683	msg = 'Expected keyword "points" but could not find it'
684	raise IOError, msg
685
686	attributes = {}
687	for key in keys:
688	if verbose: print "reading attribute '%s'" %key
689
690	attributes[key] = array(fid.variables[key])
691
692	# point_atts['attributelist'] = attributes
693
694	try:
695	geo_reference = Geo_reference(NetCDFObject=fid)
696	# point_atts['geo_reference'] = geo_reference
697	except AttributeError, e:
698	#geo_ref not compulsory
699	# point_atts['geo_reference'] = None
700	geo_reference = None
701
702	fid.close()
703
704	return pointlist, attributes, geo_reference
705
706
707	def _read_xya_file( fd, delimiter):
708	points = []
709	pointattributes = []
710	title = fd.readline()
711	att_names = clean_line(title,delimiter)
712	att_dict = {}
713	line = fd.readline()
714	numbers = clean_line(line,delimiter)
715
716	while len(numbers) > 1 and line[0] <> '#':
717	if numbers != []:
718	try:
719	x = float(numbers[0])
720	y = float(numbers[1])
721	points.append([x,y])
722	numbers.pop(0)
723	numbers.pop(0)
724	if len(att_names) != len(numbers):
725	fd.close()
726	# It might not be a problem with the title
727	#raise TitleAmountError
728	raise IOError
729	for i,num in enumerate(numbers):
730	num.strip()
731	if num != '\n' and num != '':
732	#attributes.append(float(num))
733	att_dict.setdefault(att_names[i],[]).append(float(num))
734	except ValueError:
735	raise SyntaxError
736	line = fd.readline()
737	numbers = clean_line(line,delimiter)
738
739	if line == '':
740	geo_reference = None
741	else:
742	geo_reference = Geo_reference(ASCIIFile=fd,read_title=line)
743
744
745	pointlist = array(points).astype(Float)
746	for key in att_dict.keys():
747	att_dict[key] = array(att_dict[key]).astype(Float)
748
749	return pointlist, att_dict, geo_reference
750
751	def _write_pts_file(file_name,
752	write_data_points,
753	write_attributes=None,
754	write_geo_reference=None):
755	"""
756	Write .pts NetCDF file
757
758	NOTE: Below might not be valid ask Duncan : NB 5/2006
759
760	WARNING: This function mangles the point_atts data structure
761	#F??ME: (DSG)This format has issues.
762	# There can't be an attribute called points
763	# consider format change
764	# method changed by NB not sure if above statement is correct
765
766	should create new test for this
767	legal_keys = ['pointlist', 'attributelist', 'geo_reference']
768	for key in point_atts.keys():
769	msg = 'Key %s is illegal. Valid keys are %s' %(key, legal_keys)
770	assert key in legal_keys, msg
771	"""
772	from Scientific.IO.NetCDF import NetCDFFile
773	# NetCDF file definition
774	outfile = NetCDFFile(file_name, 'w')
775
776	#Create new file
777	outfile.institution = 'Geoscience Australia'
778	outfile.description = 'NetCDF format for compact and portable storage ' +\
779	'of spatial point data'
780
781	# dimension definitions
782	shape = write_data_points.shape[0]
783	outfile.createDimension('number_of_points', shape)
784	outfile.createDimension('number_of_dimensions', 2) #This is 2d data
785
786	# variable definition
787	outfile.createVariable('points', Float, ('number_of_points',
788	'number_of_dimensions'))
789
790	#create variables
791	outfile.variables['points'][:] = write_data_points #.astype(Float32)
792
793	if write_attributes is not None:
794	for key in write_attributes.keys():
795	outfile.createVariable(key, Float, ('number_of_points',))
796	outfile.variables[key][:] = write_attributes[key] #.astype(Float32)
797
798	if write_geo_reference is not None:
799	write_geo_reference.write_NetCDF(outfile)
800
801	outfile.close()
802
803
804
805	def _write_xya_file(file_name,
806	write_data_points,
807	write_attributes=None,
808	write_geo_reference=None,
809	delimiter = ','):
810	"""
811	export a file, file_name, with the xya format
812
813	"""
814	points = write_data_points
815	pointattributes = write_attributes
816
817	fd = open(file_name,'w')
818	titlelist = ""
819	if pointattributes is not None:
820	for title in pointattributes.keys():
821	titlelist = titlelist + title + delimiter
822	titlelist = titlelist[0:-len(delimiter)] # remove the last delimiter
823	fd.write(titlelist+"\n")
824
825	#<vertex #> <x> <y> [attributes]
826	for i, vert in enumerate( points):
827
828
829	if pointattributes is not None:
830	attlist = ","
831	for att in pointattributes.keys():
832	attlist = attlist + str(pointattributes[att][i])+ delimiter
833	attlist = attlist[0:-len(delimiter)] # remove the last delimiter
834	attlist.strip()
835	else:
836	attlist = ''
837
838	fd.write(str(vert[0]) + delimiter +
839	str(vert[1]) + attlist + "\n")
840
841	if write_geo_reference is not None:
842	write_geo_reference.write_ASCII(fd)
843	fd.close()
844
845
846
847	def _point_atts2array(point_atts):
848	point_atts['pointlist'] = array(point_atts['pointlist']).astype(Float)
849
850	for key in point_atts['attributelist'].keys():
851	point_atts['attributelist'][key]= array(point_atts['attributelist'][key]).astype(Float)
852	return point_atts
853
854
855
856
857	def geospatial_data2points_dictionary(geospatial_data):
858	"""Convert geospatial data to points_dictionary
859	"""
860
861	points_dictionary = {}
862	points_dictionary['pointlist'] = geospatial_data.data_points
863
864	points_dictionary['attributelist'] = {}
865
866	for attribute_name in geospatial_data.attributes.keys():
867	val = geospatial_data.attributes[attribute_name]
868	points_dictionary['attributelist'][attribute_name] = val
869
870	points_dictionary['geo_reference'] = geospatial_data.geo_reference
871
872	return points_dictionary
873
874
875	def points_dictionary2geospatial_data(points_dictionary):
876	"""Convert points_dictionary to geospatial data object
877	"""
878
879	msg = 'Points dictionary must have key pointlist'
880	assert points_dictionary.has_key('pointlist'), msg
881
882	msg = 'Points dictionary must have key attributelist'
883	assert points_dictionary.has_key('attributelist'), msg
884
885	if points_dictionary.has_key('geo_reference'):
886	geo = points_dictionary['geo_reference']
887	else:
888	geo = None
889
890	return Geospatial_data(points_dictionary['pointlist'],
891	points_dictionary['attributelist'],
892	geo_reference = geo)
893
894	def clean_line(line,delimiter):
895	"""Remove whitespace
896	"""
897	#print ">%s" %line
898	line = line.strip()
899	#print "stripped>%s" %line
900	numbers = line.split(delimiter)
901	i = len(numbers) - 1
902	while i >= 0:
903	if numbers[i] == '':
904	numbers.pop(i)
905	else:
906	numbers[i] = numbers[i].strip()
907
908	i += -1
909	#for num in numbers:
910	# print "num>%s<" %num
911	return numbers
912
913	def ensure_absolute(points, geo_reference = None):
914	"""
915	This function inputs several formats and
916	outputs one format. - a numeric array of absolute points.
917
918	Inputed formats are;
919	points: List or numeric array of coordinate pairs [xi, eta] of
920	points or geospatial object or points file name
921
922	mesh_origin: A geo_reference object or 3-tuples consisting of
923	UTM zone, easting and northing.
924	If specified vertex coordinates are assumed to be
925	relative to their respective origins.
926	"""
927	if isinstance(points,type('')):
928	#It's a string
929	#assume it is a point file
930	points = Geospatial_data(file_name = points)
931
932	if isinstance(points,Geospatial_data):
933	points = points.get_data_points( \
934	absolute = True)
935	msg = "Use a Geospatial_data object or a mesh origin. Not both."
936	assert geo_reference == None, msg
937
938	else:
939	points = ensure_numeric(points, Float)
940	if geo_reference is None:
941	geo = None #Geo_reference()
942	else:
943	if isinstance(geo_reference, Geo_reference):
944	geo = geo_reference
945	else:
946	geo = Geo_reference(geo_reference[0],
947	geo_reference[1],
948	geo_reference[2])
949	points = geo.get_absolute(points)
950	return points
951
952
953	def ensure_geospatial(points, geo_reference = None):
954	"""
955	This function inputs several formats and
956	outputs one format. - a geospatial_data instance.
957
958	Inputed formats are;
959	points: List or numeric array of coordinate pairs [xi, eta] of
960	points or geospatial object
961
962	mesh_origin: A geo_reference object or 3-tuples consisting of
963	UTM zone, easting and northing.
964	If specified vertex coordinates are assumed to be
965	relative to their respective origins.
966	"""
967	if isinstance(points,Geospatial_data):
968	msg = "Use a Geospatial_data object or a mesh origin. Not both."
969	assert geo_reference == None, msg
970	return points
971	else:
972	points = ensure_numeric(points, Float)
973	if geo_reference is None:
974	geo = None
975	else:
976	if isinstance(geo_reference, Geo_reference):
977	geo = geo_reference
978	else:
979	geo = Geo_reference(geo_reference[0],
980	geo_reference[1],
981	geo_reference[2])
982	points = Geospatial_data(data_points=points, geo_reference=geo)
983	return points
984
985	#def file2xya(filename):
986
987	# G = Geospatial_data(filename)
988	# G.export_points_file()
989
990

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

Download in other formats: