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

Last change on this file since 3745 was 3745, checked in by ole, 17 years ago

A few minor changes to convert_from_latlon_to_UTM to make run_karratha work

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