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

Last change on this file since 3514 was 3514, checked in by duncan, 18 years ago

Hi all,
I'm doing a change in the anuga structure, moving the code to

\anuga_core\source\anuga

After you have done an svn update, the PYTHONPATH has to be changed to;
PYTHONPATH = anuga_core/source/

This is part of changes required to make installation of anuga quicker and reducing the size of our sandpits.

If any imports are broken, try fixing them. With adding anuga. to them for example. If this seems to have really broken things, email/phone me.

Cheers
Duncan

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