source: inundation/geospatial_data/geospatial_data.py @ 2711

Last change on this file since 2711 was 2711, checked in by sexton, 18 years ago

modifications for testing input

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