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

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

Added clipping method to Geospatial data.

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