Context Navigation

source: inundation/geospatial_data/geospatial_data.py @ 2698

Last change on this file since 2698 was 2698, checked in by nick, 18 years ago
updates to geospatial_data.py removed absolute from init and updated export_points_file
File size: 34.9 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
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	# absolute = True):
24
25
26	'''
27	Create instance from data points and associated attributes
28
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	'''
89
90
91	if file_name is None:
92	self.file_name = None
93	self.check_data_points(data_points)
94	self.set_attributes(attributes)
95	self.set_geo_reference(geo_reference)
96	self.set_default_attribute_name(default_attribute_name)
97	# self.set_absolute(absolute)
98
99	else:
100	if access(file_name, F_OK) == 0 :
101	msg = 'File %s does not exist or is not accessible' %file_name
102	raise msg
103	else:
104	data = {}
105	# watch for case where file name and points, attributes etc are provided!!
106	# if file name then all provided info will be removed!
107	self.file_name = file_name
108	self.new_import_points_file(file_name)
109
110	# print'data: point in init',self.data_points
111	# print'attrib: point in init',self.attributes
112	# print'geo_ref: point in init',self.geo_reference
113	self.check_data_points(self.data_points)
114	self.set_attributes(self.attributes)
115	self.set_geo_reference(self.geo_reference)
116	self.set_default_attribute_name(default_attribute_name)
117
118
119	def check_data_points(self, data_points):
120	"""Checks data points
121	"""
122
123	if data_points is None:
124	self.data_points = None
125	# FIXME: should throw an error if bad file name and no data!
126	print 'There is no data or files to read for data!!'
127	msg = 'file name %s does not exist in data set and the ' %self.file_name
128	# some assert
129
130	else:
131	self.data_points = ensure_numeric(data_points)
132
133	def set_attributes(self, attributes):
134	"""Check and assign attributes dictionary
135	"""
136
137	from types import DictType
138
139	if attributes is None:
140	self.attributes = None
141	return
142
143	if not isinstance(attributes, DictType):
144	#Convert single attribute into dictionary
145	attributes = {'attribute': attributes}
146
147	#Check input attributes
148	for key in attributes.keys():
149	try:
150	attributes[key] = ensure_numeric(attributes[key])
151	except:
152	msg = 'Attribute %s could not be converted' %key
153	msg += 'to a numeric vector'
154	raise msg
155
156	self.attributes = attributes
157
158
159	def set_geo_reference(self, geo_reference):
160
161	from coordinate_transforms.geo_reference import Geo_reference
162
163
164	if geo_reference is None:
165	self.geo_reference = Geo_reference() # Use default
166	elif isinstance(geo_reference, Geo_reference):
167	self.geo_reference = geo_reference
168	else:
169	msg = 'Argument geo_reference must be a valid Geo_reference \n'
170	msg += 'object or None.'
171	raise msg
172
173
174	def set_default_attribute_name(self, default_attribute_name):
175	self.default_attribute_name = default_attribute_name
176
177	def set_file_name(self, file_name = None):
178	if file_name is None:
179	self.file_name = None
180	else:
181	if access(file_name, F_OK) == 0 :
182	msg = 'Geospatial_data object needs have either data_points \n'
183	msg += ' or a file with data points'
184	raise msg
185	else:
186	self.file_name = file_name
187	print 'file name from set', self.file_name
188	'''
189	def set_absolute(self, absolute = True):
190	if absolute is False:
191	self.absolute = False
192	else:
193	self.absolute = True
194	'''
195	def get_geo_reference(self):
196	return self.geo_reference
197
198	def get_data_points(self, absolute = True):
199	"""Get coordinates for all data points as an Nx2 array
200
201	If absolute is True absolute UTM coordinates are returned otherwise
202	returned coordinates are relative to the internal georeference's
203	xll and yll corners.
204
205	Default: absolute is True.
206	"""
207
208	if absolute is True:
209	return self.geo_reference.get_absolute(self.data_points)
210	else:
211	return self.data_points
212
213
214	def get_attributes(self, attribute_name = None):
215	"""Return values for one named attribute.
216
217	If attribute_name is None, default_attribute_name is used
218	"""
219
220	if attribute_name is None:
221	if self.default_attribute_name is not None:
222	attribute_name = self.default_attribute_name
223	else:
224	attribute_name = self.attributes.keys()[0]
225	# above line takes the first one from keys
226
227
228	msg = 'Attribute name %s does not exist in data set' %attribute_name
229	assert self.attributes.has_key(attribute_name), msg
230
231	return self.attributes[attribute_name]
232
233	def get_all_attributes(self):
234	"""
235	Return values for all attributes.
236	"""
237
238	return self.attributes
239
240	def __add__(self, other):
241	'''
242	returns the addition of 2 geospatical objects,
243	objects are concatencated to the end of each other
244
245	NOTE: doesn't add if objects contain different
246	attributes
247
248	Always return relative points!
249	'''
250
251
252	# find objects zone and checks if the same
253	geo_ref1 = self.get_geo_reference()
254	zone1 = geo_ref1.get_zone()
255
256	geo_ref2 = other.get_geo_reference()
257	zone2 = geo_ref2.get_zone()
258
259	geo_ref1.reconcile_zones(geo_ref2)
260
261
262	# sets xll and yll as the smallest from self and other
263	# FIXME (Duncan and Ole): use lower left corner derived from absolute coordinates
264	if self.geo_reference.xllcorner <= other.geo_reference.xllcorner:
265	xll = self.geo_reference.xllcorner
266	else:
267	xll = other.geo_reference.xllcorner
268
269	if self.geo_reference.yllcorner <= other.geo_reference.yllcorner:
270	yll = self.geo_reference.yllcorner
271	else:
272	yll = other.geo_reference.yllcorner
273	new_geo_ref = Geo_reference(geo_ref1.get_zone(), xll, yll)
274
275	xll = yll = 0.
276
277	relative_points1 = self.get_data_points(absolute = False)
278	relative_points2 = other.get_data_points(absolute = False)
279
280
281	new_relative_points1 = new_geo_ref.change_points_geo_ref(relative_points1, geo_ref1)
282	new_relative_points2 = new_geo_ref.change_points_geo_ref(relative_points2, geo_ref2)
283
284	#Now both point sets are relative to new_geo_ref and zones have been reconciled
285
286
287	# Concatenate points
288	new_points = concatenate((new_relative_points1,
289	new_relative_points2),
290	axis = 0)
291	# print 'new points:', new_points
292
293
294	# Concatenate attributes
295	new_attributes = {}
296	for x in self.attributes.keys():
297	if other.attributes.has_key(x):
298
299	attrib1 = self.attributes[x]
300	attrib2 = other.attributes[x]
301	new_attributes[x] = concatenate((attrib1, attrib2))
302
303	else:
304	msg = 'Both geospatial_data objects must have the same \n'
305	msg += 'attributes to allow addition.'
306	raise msg
307
308
309	# Instantiate new data object and return
310	return Geospatial_data(new_points,
311	new_attributes,
312	new_geo_ref)
313
314
315	'''
316	def xxx__add__(self, other):
317	"""
318	returns the addition of 2 geospatical objects,
319	objects are concatencated to the end of each other
320
321	NOTE: doesn't add if objects contain different
322	attributes
323	"""
324
325	# sets xll and yll as the smallest from self and other
326	if self.geo_reference.xllcorner <= other.geo_reference.xllcorner:
327	xll = self.geo_reference.xllcorner
328	else:
329	xll = other.geo_reference.xllcorner
330
331	if self.geo_reference.yllcorner <= other.geo_reference.yllcorner:
332	yll = self.geo_reference.yllcorner
333	else:
334	yll = other.geo_reference.yllcorner
335
336	# find objects zone and checks if the same
337	geo_ref1 = self.get_geo_reference()
338	zone1 = geo_ref1.get_zone()
339
340	geo_ref2 = other.get_geo_reference()
341	zone2 = geo_ref2.get_zone()
342
343	if zone1 == zone2:
344	a_absolute_points = self.get_data_points(absolute=True)
345	b_absolute_points = other.get_data_points(absolute=True)
346
347	# subtract xll and yll from self's and other's
348	# absolute data points and concatenate
349	c_points = concatenate((a_absolute_points-[xll, yll],
350	b_absolute_points-[xll, yll]),
351	axis = 0)
352
353	new_attributes = {}
354	for x in self.attributes.keys():
355	if other.attributes.has_key(x):
356
357	a_attrib = self.attributes[x]
358	b_attrib = other.attributes[x]
359	new_attributes[x] = concatenate((a_attrib, b_attrib))
360
361	else:
362	msg = 'Both geospatial_data objects must have the same \n'
363	msg += 'attributes to allow addition.'
364	raise msg
365
366	geo_ref = Geo_reference(zone1, xll, yll)
367	return Geospatial_data(c_points, new_attributes, geo_ref)
368
369	else:
370	msg = 'Both geospatial_data objects must be in the same \
371	ZONE to allow addition. I got zone %d and %d' %(zone1, zone2)
372	raise msg
373
374	'''
375
376	###
377	# IMPORT/EXPORT POINTS FILES
378	###
379	def new_import_points_file(self, ofile, delimiter = None, verbose = False):
380	""" load an .xya or .pts file
381	Note: will throw an IOError if it can't load the file.
382	Catch these!
383	"""
384
385	attributes = {}
386	if ofile[-4:]== ".xya":
387	try:
388	if delimiter == None:
389	try:
390	fd = open(ofile)
391	# all_data = _read_xya_file(fd, ',')
392	data_points, attributes, geo_reference = _new_read_xya_file(fd, ',')
393	except SyntaxError:
394	fd.close()
395	fd = open(ofile)
396	# all_data = _read_xya_file(fd, ' ')
397	data_points, attributes, geo_reference = _new_read_xya_file(fd, ' ')
398	else:
399	fd = open(ofile)
400	# all_data = _read_xya_file(fd, delimiter)
401	data_points, attributes, geo_reference = _new_read_xya_file(fd, delimiter)
402	fd.close()
403	except (IndexError,ValueError,SyntaxError):
404	fd.close()
405	msg = 'Could not open file %s ' %ofile
406	raise IOError, msg
407	except IOError, e:
408	# Catch this to add an error message
409	msg = 'Could not open file or incorrect file format %s:%s' %(ofile, e)
410	raise IOError, msg
411
412	elif ofile[-4:]== ".pts":
413	try:
414	# print 'hi from import_points_file'
415	# all_data = _read_pts_file(ofile, verbose)
416	data_points, attributes, geo_reference = _new_read_pts_file(ofile, verbose)
417	# print 'hi1 from import_points_file', all_data
418	except IOError, e:
419	msg = 'Could not open file %s ' %ofile
420	raise IOError, msg
421	else:
422	msg = 'Extension %s is unknown' %ofile[-4:]
423	raise IOError, msg
424
425	# print'in import data_points', data_points
426	# print'in import attributes', attributes
427	# print'in import data_points', geo_reference
428	self.data_points = data_points
429	self.attributes = attributes
430	self.geo_reference = geo_reference
431
432	# return all_data
433
434	def export_points_file(self, ofile, absolute=False):
435
436	"""
437	write a points file, ofile, as a text (.xya) or binary (.pts) file
438	ofile is the file name, including the extension
439	The point_dict is defined at the top of this file.
440	"""
441
442	#this was done for all keys in the mesh file.
443	#if not mesh_dict.has_key('points'):
444	# mesh_dict['points'] = []
445
446	if (ofile[-4:] == ".xya"):
447	if absolute is True:
448	_write_xya_file(ofile, self.get_data_points(absolute),
449	self.get_all_attributes())
450	else:
451	_write_xya_file(ofile, self.get_data_points(absolute),
452	self.get_all_attributes(),
453	self.get_geo_reference())
454
455	elif (ofile[-4:] == ".pts"):
456	if absolute is True:
457	_write_pts_file(ofile, self.get_data_points(absolute),
458	self.get_all_attributes())
459	else:
460	_write_pts_file(ofile, self.get_data_points(absolute),
461	self.get_all_attributes(),
462	self.get_geo_reference())
463	else:
464	msg = 'Unknown file type %s ' %ofile
465	raise IOError, msg
466
467
468	def import_points_file( ofile, delimiter = None, verbose = False):
469	""" load an .xya or .pts file
470	Note: will throw an IOError if it can't load the file.
471	Catch these!
472	"""
473
474	all_data = {}
475	if ofile[-4:]== ".xya":
476	try:
477	if delimiter == None:
478	try:
479	fd = open(ofile)
480	all_data = _read_xya_file(fd, ',')
481	except SyntaxError:
482	fd.close()
483	fd = open(ofile)
484	all_data = _read_xya_file(fd, ' ')
485	else:
486	fd = open(ofile)
487	all_data = _read_xya_file(fd, delimiter)
488	fd.close()
489	except (IndexError,ValueError,SyntaxError):
490	fd.close()
491	msg = 'Could not open file %s ' %ofile
492	raise IOError, msg
493	except IOError, e:
494	# Catch this to add an error message
495	msg = 'Could not open file or incorrect file format %s:%s' %(ofile, e)
496	raise IOError, msg
497
498	elif ofile[-4:]== ".pts":
499	try:
500	# print 'hi from import_points_file'
501	all_data = _read_pts_file(ofile, verbose)
502	# print 'hi1 from import_points_file', all_data
503	except IOError, e:
504	msg = 'Could not open file %s ' %ofile
505	raise IOError, msg
506	else:
507	msg = 'Extension %s is unknown' %ofile[-4:]
508	raise IOError, msg
509
510	return all_data
511	'''
512	def export_points_file(ofile, point_dict):
513	"""
514	write a points file, ofile, as a text (.xya) or binary (.pts) file
515
516	ofile is the file name, including the extension
517
518	The point_dict is defined at the top of this file.
519	"""
520	#this was done for all keys in the mesh file.
521	#if not mesh_dict.has_key('points'):
522	# mesh_dict['points'] = []
523	if (ofile[-4:] == ".xya"):
524	_write_xya_file(ofile, point_dict)
525	elif (ofile[-4:] == ".pts"):
526	_write_pts_file(ofile, point_dict)
527	else:
528	msg = 'Unknown file type %s ' %ofile
529	raise IOError, msg
530	'''
531
532
533	def _new_read_pts_file(file_name, verbose = False):
534	"""Read .pts NetCDF file
535
536	Return a dic of array of points, and dic of array of attribute
537	eg
538	dic['points'] = [[1.0,2.0],[3.0,5.0]]
539	dic['attributelist']['elevation'] = [[7.0,5.0]
540	"""
541	#FIXME: (DSG) This format has issues.
542	# There can't be an attribute called points
543	# consider format change
544
545	# print 'hi for read points file'
546	from Scientific.IO.NetCDF import NetCDFFile
547
548	if verbose: print 'Reading ', file_name
549
550
551	# see if the file is there. Throw a QUIET IO error if it isn't
552	fd = open(file_name,'r')
553	fd.close()
554
555	#throws prints to screen if file not present
556	fid = NetCDFFile(file_name, 'r')
557
558	# point_atts = {}
559	# Get the variables
560	# point_atts['pointlist'] = array(fid.variables['points'])
561	pointlist = array(fid.variables['points'])
562	keys = fid.variables.keys()
563	if verbose: print 'Got %d variables: %s' %(len(keys), keys)
564	try:
565	keys.remove('points')
566	except IOError, e:
567	fid.close()
568	msg = 'Expected keyword "points" but could not find it'
569	raise IOError, msg
570
571	attributes = {}
572	for key in keys:
573	if verbose: print "reading attribute '%s'" %key
574
575	attributes[key] = array(fid.variables[key])
576
577	# point_atts['attributelist'] = attributes
578
579	try:
580	geo_reference = Geo_reference(NetCDFObject=fid)
581	# point_atts['geo_reference'] = geo_reference
582	except AttributeError, e:
583	#geo_ref not compulsory
584	# point_atts['geo_reference'] = None
585	geo_reference = None
586
587	fid.close()
588
589	return pointlist, attributes, geo_reference
590
591
592	def _new_read_xya_file( fd, delimiter):
593	# print 'hello from read xya data'
594	points = []
595	pointattributes = []
596	title = fd.readline()
597	att_names = clean_line(title,delimiter)
598
599	att_dict = {}
600	line = fd.readline()
601	numbers = clean_line(line,delimiter)
602	#print 'read xya numbers', numbers
603	while len(numbers) > 1:
604	if numbers != []:
605	try:
606	x = float(numbers[0])
607	y = float(numbers[1])
608	points.append([x,y])
609	numbers.pop(0)
610	numbers.pop(0)
611	#attributes = []
612	#print "att_names",att_names
613	#print "numbers",numbers
614	if len(att_names) != len(numbers):
615	fd.close()
616	# It might not be a problem with the title
617	#raise TitleAmountError
618	raise IOError
619	for i,num in enumerate(numbers):
620	num.strip()
621	if num != '\n' and num != '':
622	#attributes.append(float(num))
623	att_dict.setdefault(att_names[i],[]).append(float(num))
624
625	except ValueError:
626	raise SyntaxError
627	line = fd.readline()
628	numbers = clean_line(line,delimiter)
629
630	if line == '':
631	# end of file
632	geo_reference = None
633	else:
634	geo_reference = Geo_reference(ASCIIFile=fd,read_title=line)
635
636	# xya_dict = {}
637	# xya_dict['pointlist'] = array(points).astype(Float)
638	pointlist = array(points).astype(Float)
639
640	for key in att_dict.keys():
641	att_dict[key] = array(att_dict[key]).astype(Float)
642	# xya_dict['attributelist'] = att_dict
643
644	# xya_dict['geo_reference'] = geo_reference
645	#print "xya_dict",xya_dict
646	return pointlist, att_dict, geo_reference
647
648
649	def _read_pts_file(file_name, verbose = False):
650	"""Read .pts NetCDF file
651
652	Return a dic of array of points, and dic of array of attribute
653	eg
654	dic['points'] = [[1.0,2.0],[3.0,5.0]]
655	dic['attributelist']['elevation'] = [[7.0,5.0]
656	"""
657	#FIXME: (DSG) This format has issues.
658	# There can't be an attribute called points
659	# consider format change
660
661	# print 'hi for read points file'
662	from Scientific.IO.NetCDF import NetCDFFile
663
664	if verbose: print 'Reading ', file_name
665
666
667	# see if the file is there. Throw a QUIET IO error if it isn't
668	fd = open(file_name,'r')
669	fd.close()
670
671	#throws prints to screen if file not present
672	fid = NetCDFFile(file_name, 'r')
673
674	point_atts = {}
675	# Get the variables
676	point_atts['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
701	fid.close()
702	return point_atts
703
704
705	def _read_xya_file( fd, delimiter):
706	# print 'hello from read xya data'
707	points = []
708	pointattributes = []
709	title = fd.readline()
710	att_names = clean_line(title,delimiter)
711
712	att_dict = {}
713	line = fd.readline()
714	numbers = clean_line(line,delimiter)
715	#print 'read xya numbers', numbers
716	while len(numbers) > 1:
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	#attributes = []
725	#print "att_names",att_names
726	#print "numbers",numbers
727	if len(att_names) != len(numbers):
728	fd.close()
729	# It might not be a problem with the title
730	#raise TitleAmountError
731	raise IOError
732	for i,num in enumerate(numbers):
733	num.strip()
734	if num != '\n' and num != '':
735	#attributes.append(float(num))
736	att_dict.setdefault(att_names[i],[]).append(float(num))
737
738	except ValueError:
739	raise SyntaxError
740	line = fd.readline()
741	numbers = clean_line(line,delimiter)
742
743	if line == '':
744	# end of file
745	geo_reference = None
746	else:
747	geo_reference = Geo_reference(ASCIIFile=fd,read_title=line)
748
749	xya_dict = {}
750	xya_dict['pointlist'] = array(points).astype(Float)
751
752	for key in att_dict.keys():
753	att_dict[key] = array(att_dict[key]).astype(Float)
754	xya_dict['attributelist'] = att_dict
755
756	xya_dict['geo_reference'] = geo_reference
757	#print "xya_dict",xya_dict
758	return xya_dict
759
760
761	def _write_pts_file(file_name, write_data_points,
762	write_attributes,
763	write_geo_reference = None):
764	"""Write .pts NetCDF file
765
766	WARNING: This function mangles the point_atts data structure
767	"""
768	#FIXME: (DSG)This format has issues.
769	# There can't be an attribute called points
770	# consider format change
771	# method changed by NB not sure if above statement is correct
772
773	''' should create new test for this
774	legal_keys = ['pointlist', 'attributelist', 'geo_reference']
775	for key in point_atts.keys():
776	msg = 'Key %s is illegal. Valid keys are %s' %(key, legal_keys)
777	assert key in legal_keys, msg
778	'''
779	from Scientific.IO.NetCDF import NetCDFFile
780	# NetCDF file definition
781	outfile = NetCDFFile(file_name, 'w')
782
783	#Create new file
784	outfile.institution = 'Geoscience Australia'
785	outfile.description = 'NetCDF format for compact and portable storage ' +\
786	'of spatial point data'
787
788	# dimension definitions
789	shape = write_data_points.shape[0]
790	outfile.createDimension('number_of_points', shape)
791	outfile.createDimension('number_of_dimensions', 2) #This is 2d data
792
793	# variable definition
794	outfile.createVariable('points', Float, ('number_of_points',
795	'number_of_dimensions'))
796
797	#create variables
798	outfile.variables['points'][:] = write_data_points #.astype(Float32)
799
800	for key in write_attributes.keys():
801	outfile.createVariable(key, Float, ('number_of_points',))
802	outfile.variables[key][:] = write_attributes[key] #.astype(Float32)
803
804	if write_geo_reference is not None:
805	write_geo_reference.write_NetCDF(outfile)
806
807	outfile.close()
808
809
810
811	def _write_xya_file( file_name, write_data_points,
812	write_attributes,
813	write_geo_reference = None,
814	delimiter = ','):
815	"""
816	export a file, ofile, with the xya format
817
818	"""
819	points = write_data_points
820	pointattributes = write_attributes
821
822	fd = open(file_name,'w')
823	titlelist = ""
824	for title in pointattributes.keys():
825	titlelist = titlelist + title + delimiter
826	titlelist = titlelist[0:-len(delimiter)] # remove the last delimiter
827	fd.write(titlelist+"\n")
828	#<vertex #> <x> <y> [attributes]
829	for i,vert in enumerate( points):
830
831	attlist = ","
832	for att in pointattributes.keys():
833	attlist = attlist + str(pointattributes[att][i])+ delimiter
834	attlist = attlist[0:-len(delimiter)] # remove the last delimiter
835	attlist.strip()
836	fd.write( str(vert[0]) + delimiter
837	+ str(vert[1])
838	+ attlist + "\n")
839	'''
840	# geo_reference info
841	if xya_dict.has_key('geo_reference') and \
842	not xya_dict['geo_reference'] is None:
843	xya_dict['geo_reference'].write_ASCII(fd)
844	'''
845	if write_geo_reference is not None:
846	write_geo_reference.write_ASCII(fd)
847	fd.close()
848	'''
849	def _write_pts_file(file_name, point_atts):
850	"""Write .pts NetCDF file
851
852	WARNING: This function mangles the point_atts data structure
853	"""
854	#FIXME: (DSG)This format has issues.
855	# There can't be an attribute called points
856	# consider format change
857
858
859	legal_keys = ['pointlist', 'attributelist', 'geo_reference']
860	for key in point_atts.keys():
861	msg = 'Key %s is illegal. Valid keys are %s' %(key, legal_keys)
862	assert key in legal_keys, msg
863
864	from Scientific.IO.NetCDF import NetCDFFile
865	_point_atts2array(point_atts)
866	# NetCDF file definition
867	outfile = NetCDFFile(file_name, 'w')
868
869	#Create new file
870	outfile.institution = 'Geoscience Australia'
871	outfile.description = 'NetCDF format for compact and portable storage ' +\
872	'of spatial point data'
873
874	# dimension definitions
875	shape = point_atts['pointlist'].shape[0]
876	outfile.createDimension('number_of_points', shape)
877	outfile.createDimension('number_of_dimensions', 2) #This is 2d data
878
879	# variable definition
880	outfile.createVariable('points', Float, ('number_of_points',
881	'number_of_dimensions'))
882
883	#create variables
884	outfile.variables['points'][:] = point_atts['pointlist'] #.astype(Float32)
885	for key in point_atts['attributelist'].keys():
886	outfile.createVariable(key, Float, ('number_of_points',))
887	outfile.variables[key][:] = point_atts['attributelist'][key] #.astype(Float32)
888
889	if point_atts.has_key('geo_reference') and not point_atts['geo_reference'] == None:
890	point_atts['geo_reference'].write_NetCDF(outfile)
891
892	outfile.close()
893
894
895
896	def _write_xya_file( file_name, xya_dict, delimiter = ','):
897	"""
898	export a file, ofile, with the xya format
899
900	"""
901	points = xya_dict['pointlist']
902	pointattributes = xya_dict['attributelist']
903
904	fd = open(file_name,'w')
905
906	titlelist = ""
907	for title in pointattributes.keys():
908	titlelist = titlelist + title + delimiter
909	titlelist = titlelist[0:-len(delimiter)] # remove the last delimiter
910	fd.write(titlelist+"\n")
911	#<vertex #> <x> <y> [attributes]
912	for i,vert in enumerate( points):
913
914	attlist = ","
915	for att in pointattributes.keys():
916	attlist = attlist + str(pointattributes[att][i])+ delimiter
917	attlist = attlist[0:-len(delimiter)] # remove the last delimiter
918	attlist.strip()
919	fd.write( str(vert[0]) + delimiter
920	+ str(vert[1])
921	+ attlist + "\n")
922
923	# geo_reference info
924	if xya_dict.has_key('geo_reference') and \
925	not xya_dict['geo_reference'] is None:
926	xya_dict['geo_reference'].write_ASCII(fd)
927	fd.close()
928	'''
929	def _point_atts2array(point_atts):
930	point_atts['pointlist'] = array(point_atts['pointlist']).astype(Float)
931
932	for key in point_atts['attributelist'].keys():
933	point_atts['attributelist'][key]= array(point_atts['attributelist'][key]).astype(Float)
934	return point_atts
935
936
937
938
939	def geospatial_data2points_dictionary(geospatial_data):
940	"""Convert geospatial data to points_dictionary
941	"""
942
943	points_dictionary = {}
944	points_dictionary['pointlist'] = geospatial_data.data_points
945
946	points_dictionary['attributelist'] = {}
947
948	for attribute_name in geospatial_data.attributes.keys():
949	val = geospatial_data.attributes[attribute_name]
950	points_dictionary['attributelist'][attribute_name] = val
951
952	points_dictionary['geo_reference'] = geospatial_data.geo_reference
953
954	return points_dictionary
955
956
957	def points_dictionary2geospatial_data(points_dictionary):
958	"""Convert points_dictionary to geospatial data object
959	"""
960
961	msg = 'Points dictionary must have key pointlist'
962	assert points_dictionary.has_key('pointlist'), msg
963
964	msg = 'Points dictionary must have key attributelist'
965	assert points_dictionary.has_key('attributelist'), msg
966
967	if points_dictionary.has_key('geo_reference'):
968	geo = points_dictionary['geo_reference']
969	else:
970	geo = None
971
972	return Geospatial_data(points_dictionary['pointlist'],
973	points_dictionary['attributelist'],
974	geo_reference = geo)
975
976	def clean_line(line,delimiter):
977	"""Remove whitespace
978	"""
979	#print ">%s" %line
980	line = line.strip()
981	#print "stripped>%s" %line
982	numbers = line.split(delimiter)
983	i = len(numbers) - 1
984	while i >= 0:
985	if numbers[i] == '':
986	numbers.pop(i)
987	else:
988	numbers[i] = numbers[i].strip()
989
990	i += -1
991	#for num in numbers:
992	# print "num>%s<" %num
993	return numbers
994
995	def add_points_files(add_file1, add_file2, results_file):
996	''' adds the points and attruibutes of 2 pts or xya files and
997	writes it to a pts file
998
999	NOTE will add the function to check and remove points from one set
1000	that are shared. This will require some work and maybe __subtract__ function
1001	'''
1002
1003	G1 = Geospatial_data(file_name = add_file1)
1004	G2 = Geospatial_data(file_name = add_file2)
1005	new_add_file2 = add_file2[:-4] + '.pts'
1006
1007	G = G1 + G2
1008
1009	#FIXME remove dependance on points to dict in export only!
1010	# G_points_dict = geospatial_data2points_dictionary(G)
1011	# export_points_file(results_file, G_points_dict)
1012
1013	# G_points_dict = geospatial_data2points_dictionary(G)
1014
1015	G.export_points_file(results_file)
1016
1017
1018
1019

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

Download in other formats: