Context Navigation

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

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

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

Download in other formats: