Context Navigation

source: inundation/geospatial_data/geospatial_data.py @ 3055

Last change on this file since 3055 was 3051, checked in by duncan, 19 years ago
checking in changes half way thru, so Ole can help
File size: 24.3 KB

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

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

Download in other formats: