Context Navigation

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

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

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

Download in other formats: