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

Last change on this file since 4641 was 4641, checked in by nick, 18 years ago

updated geospatial_data.export_data_points_file so it will write data "as_lat_long" in the northern hemisphere to csv and txt files.

Plus change parameter name to be more logical

and corrected a verbose=True in test_geospatial_data.py, Dunc and i didn't know why it was there...

File size: 55.3 KB
Line 
1"""Class Geospatial_data - Manipulation of locations on the planet and
2associated attributes.
3
4"""
5from sys import maxint
6from os import access, F_OK, R_OK
7from types import DictType
8from warnings import warn
9from string import lower
10from Numeric import concatenate, array, Float, shape, reshape, ravel, take, \
11                        size, shape
12#from Array import tolist
13from RandomArray import randint
14from copy import deepcopy
15
16#from MA import tolist
17
18from Scientific.IO.NetCDF import NetCDFFile   
19from anuga.coordinate_transforms.lat_long_UTM_conversion import UTMtoLL   
20from anuga.utilities.numerical_tools import ensure_numeric
21from anuga.coordinate_transforms.geo_reference import Geo_reference, \
22     TitleError, DEFAULT_ZONE, ensure_geo_reference, write_NetCDF_georeference
23from anuga.coordinate_transforms.redfearn import convert_from_latlon_to_utm
24from anuga.utilities.anuga_exceptions import ANUGAError
25from anuga.config import points_file_block_line_size as MAX_READ_LINES   
26
27DEFAULT_ATTRIBUTE = 'elevation'
28
29class Geospatial_data:
30
31    def __init__(self,
32                 data_points=None, # this can also be a points file name
33                 attributes=None,
34                 geo_reference=None,
35                 default_attribute_name=None,
36                 file_name=None,
37                 delimiter=None,
38                 latitudes=None,
39                 longitudes=None,
40                 points_are_lats_longs=False,
41                 max_read_lines=None,
42                 load_file_now=True,
43                 verbose=False):
44
45       
46        """
47        Create instance from data points and associated attributes
48
49        data_points: x,y coordinates in meters. Type must be either a
50        sequence of 2-tuples or an Mx2 Numeric array of floats.  A file name
51        with extension .txt, .cvs or .pts can also be passed in here.
52
53        attributes: Associated values for each data point. The type
54        must be either a list or an array of length M or a dictionary
55        of lists (or arrays) of length M. In the latter case the keys
56        in the dictionary represent the attribute names, in the former
57        the attribute will get the default name "elevation".
58       
59        geo_reference: Object representing the origin of the data
60        points. It contains UTM zone, easting and northing and data
61        points are assumed to be relative to this origin.
62        If geo_reference is None, the default geo ref object is used.
63
64        default_attribute_name: Name of default attribute to be used with
65        get_attribute_values. The idea is that the dataset can be
66        equipped with information about which attribute to return.
67        If None, the default is the "first"
68
69        latitudes, longitudes: Vectors of latitudes and longitudes,
70        used to specify location instead of points.
71       
72        points_are_lats_longs: Set this as true if the points are actually
73        lats and longs, not UTM
74
75        max_read_lines: The number of rows read into memory when using
76        blocking to read a file.
77
78        load_file_now:  If true the file is automatically loaded
79        into the geospatial instance. Used when blocking.
80       
81        file_name: Name of input netCDF file or .txt file. netCDF file must
82        have dimensions "points" etc.
83        .txt file is a comma seperated file with x, y and attribute
84        data.
85       
86        The first line has the titles of the columns.  The first two
87        column titles are checked to see if they start with lat or
88        long (not case sensitive).  If so the data is assumed to be
89        latitude and longitude, in decimal format and converted to
90        UTM.  Otherwise the first two columns are assumed to be the x
91        and y, and the title names acually used are ignored.
92 
93       
94        The format for a .txt file is:
95            1st line:     [column names]
96            other lines:  x y [attributes]
97
98            for example:
99            x, y, elevation, friction
100            0.6, 0.7, 4.9, 0.3
101            1.9, 2.8, 5, 0.3
102            2.7, 2.4, 5.2, 0.3
103
104        The first two columns are always  assumed to be x, y
105        coordinates.
106     
107        An issue with the xya format is that the attribute column order
108        is not be controlled.  The info is stored in a dictionary and it's
109        written in an order dependent on the hash order
110       
111        The format for a Points dictionary is:
112
113          ['pointlist'] a 2 column array describing points. 1st column x,
114          2nd column y.
115          ['attributelist'], a dictionary of 1D arrays, representing
116          attribute values at the point.  The dictionary key is the attribute
117          header.
118          ['geo_reference'] a Geo_refernece object. Use if the point
119          information is relative. This is optional.
120            eg
121            dic['pointlist'] = [[1.0,2.0],[3.0,5.0]]
122            dic['attributelist']['elevation'] = [[7.0,5.0]
123               
124        delimiter: is the file delimiter that will be used when
125            importing a .xya file, which is being phased out.
126           
127        verbose:
128
129         
130        """
131
132        if isinstance(data_points, basestring):
133            # assume data point is really a file name
134            file_name = data_points
135
136        self.set_verbose(verbose)
137        self.geo_reference=None #create the attribute
138        self.file_name = file_name
139        self.max_read_lines = max_read_lines
140
141        if delimiter is not None:
142            msg = 'Specifying delimiters will be removed.'
143            msg = 'Text file format is moving to comma seperated .txt files.'
144            warn(msg, DeprecationWarning) 
145        if file_name is None:
146            if delimiter is not None:
147                msg = 'No file specified yet a delimiter is provided!'
148                raise ValueError, msg
149           
150            if latitudes is not None or longitudes is not None or \
151                   points_are_lats_longs:
152                data_points, geo_reference =  \
153                             _set_using_lat_long(latitudes=latitudes,
154                                  longitudes=longitudes,
155                                  geo_reference=geo_reference,
156                                  data_points=data_points,
157                                  points_are_lats_longs=points_are_lats_longs)
158            self.check_data_points(data_points)
159            self.set_attributes(attributes)
160            self.set_geo_reference(geo_reference)
161            self.set_default_attribute_name(default_attribute_name)
162
163        elif load_file_now is True:
164            # watch for case where file name and points,
165            # attributes etc are provided!!
166            # if file name then all provided info will be removed!
167
168            if verbose is True:
169                if file_name is not None:
170                    print 'Loading Geospatial data from file: %s' %file_name
171           
172            self.import_points_file(file_name, delimiter, verbose)
173               
174            self.check_data_points(self.data_points)
175            self.set_attributes(self.attributes) 
176            self.set_geo_reference(self.geo_reference)
177            self.set_default_attribute_name(default_attribute_name)
178
179
180        #Why?   
181        #assert self.attributes is None or isinstance(self.attributes, DictType)
182        #This is a hassle when blocking, so I've removed it.
183
184
185        if verbose is True:
186            if file_name is not None:
187                print 'Geospatial data created from file: %s' %file_name
188                if load_file_now is False:
189                    print 'Data will be loaded blockwise on demand'
190
191                    if file_name.endswith('csv') or file_name.endswith('txt'):
192                        print 'ASCII formats are not that great for '
193                        print 'blockwise reading. Consider storing this'
194                        print 'data as a pts NetCDF format'
195
196    def __len__(self):
197        return len(self.data_points)
198
199    def __repr__(self):
200        return str(self.get_data_points(absolute=True))
201   
202   
203    def check_data_points(self, data_points):
204        """Checks data points
205        """
206       
207        if data_points is None:
208            self.data_points = None
209            msg = 'There is no data or file provided!'
210            raise ValueError, msg
211           
212        else:
213            self.data_points = ensure_numeric(data_points)
214            #print "self.data_points.shape",self.data_points.shape
215            if not (0,) == self.data_points.shape:
216                assert len(self.data_points.shape) == 2
217                assert self.data_points.shape[1] == 2
218
219    def set_attributes(self, attributes):
220        """Check and assign attributes dictionary
221        """
222       
223        if attributes is None:
224            self.attributes = None
225            return
226       
227        if not isinstance(attributes, DictType):
228            #Convert single attribute into dictionary
229            attributes = {DEFAULT_ATTRIBUTE: attributes}
230
231        #Check input attributes   
232        for key in attributes.keys():
233            try:
234                attributes[key] = ensure_numeric(attributes[key])
235            except:
236                msg = 'Attribute %s could not be converted' %key
237                msg += 'to a numeric vector'
238                raise msg
239
240        self.attributes = attributes   
241
242
243    def set_geo_reference(self, geo_reference):
244        """
245        Set's the georeference of geospatial.
246        It can also be used to change the georeference
247        """
248        from anuga.coordinate_transforms.geo_reference import Geo_reference
249
250        if geo_reference is None:
251            geo_reference = Geo_reference() # Use default
252        geo_reference = ensure_geo_reference(geo_reference)
253        if not isinstance(geo_reference, Geo_reference):
254            msg = 'Argument geo_reference must be a valid Geo_reference \n'
255            msg += 'object or None.'
256            raise msg
257
258        # if a geo ref already exists, change the point data to
259        # represent the new geo-ref
260        if  self.geo_reference is not None:
261            #FIXME: Maybe put out a warning here...
262            self.data_points = self.get_data_points \
263                               (geo_reference=geo_reference)
264           
265        self.geo_reference = geo_reference
266
267
268    def set_default_attribute_name(self, default_attribute_name):
269        self.default_attribute_name = default_attribute_name
270
271    def set_verbose(self, verbose=False):
272        if verbose in [False, True]:
273            self.verbose = verbose
274        else:
275            msg = 'Illegal value: %s' %str(verbose)
276            raise Exception, msg
277
278    def clip(self, polygon, closed=True):
279        """Clip geospatial data by a polygon
280
281        Input
282          polygon - Either a list of points, an Nx2 array or
283                    a Geospatial data object.
284          closed - (optional) determine whether points on boundary should be
285          regarded as belonging to the polygon (closed = True)
286          or not (closed = False). Default is True.
287           
288        Output
289          New geospatial data object representing points inside
290          specified polygon.
291       
292        """
293
294        from anuga.utilities.polygon import inside_polygon
295
296        if isinstance(polygon, Geospatial_data):
297            # Polygon is an object - extract points
298            polygon = polygon.get_data_points()
299
300        points = self.get_data_points()   
301        inside_indices = inside_polygon(points, polygon, closed)
302
303        clipped_G = self.get_sample(inside_indices)
304#        clipped_points = take(points, inside_indices)
305
306        # Clip all attributes
307#        attributes = self.get_all_attributes()
308
309#        clipped_attributes = {}
310#        if attributes is not None:
311#            for key, att in attributes.items():
312#                clipped_attributes[key] = take(att, inside_indices)
313
314#        return Geospatial_data(clipped_points, clipped_attributes)
315        return clipped_G
316       
317       
318    def clip_outside(self, polygon, closed=True):
319        """Clip geospatial date by a polygon, keeping data OUTSIDE of polygon
320
321        Input
322          polygon - Either a list of points, an Nx2 array or
323                    a Geospatial data object.
324          closed - (optional) determine whether points on boundary should be
325          regarded as belonging to the polygon (closed = True)
326          or not (closed = False). Default is True.
327         
328        Output
329          Geospatial data object representing point OUTSIDE specified polygon
330       
331        """
332
333        from anuga.utilities.polygon import outside_polygon
334
335        if isinstance(polygon, Geospatial_data):
336            # Polygon is an object - extract points
337            polygon = polygon.get_data_points()
338
339        points = self.get_data_points()   
340        outside_indices = outside_polygon(points, polygon, closed)
341
342        clipped_G = self.get_sample(outside_indices)
343
344#        clipped_points = take(points, outside_indices)
345
346        # Clip all attributes
347#        attributes = self.get_all_attributes()
348
349#        clipped_attributes = {}
350#        if attributes is not None:
351#            for key, att in attributes.items():
352#                clipped_attributes[key] = take(att, outside_indices)
353
354#        return Geospatial_data(clipped_points, clipped_attributes)
355        return clipped_G
356
357   
358    def get_geo_reference(self):
359        return self.geo_reference
360       
361    def get_data_points(self, absolute=True, geo_reference=None,
362                        as_lat_long=False,isSouthHemisphere=True):
363        """Get coordinates for all data points as an Nx2 array
364
365        If absolute is False returned coordinates are relative to the
366        internal georeference's xll and yll corners, otherwise
367        absolute UTM coordinates are returned.
368
369        If a geo_reference is passed the points are returned relative
370        to that geo_reference.
371       
372        isSH (isSouthHemisphere) is only used when getting data
373        points "as_lat_long" is True and if FALSE will return lats and
374        longs valid for the Northern Hemisphere.
375
376        Default: absolute is True.
377        """
378        if as_lat_long is True:
379            msg = "Points need a zone to be converted into lats and longs"
380            assert self.geo_reference is not None, msg
381            zone = self.geo_reference.get_zone()
382            assert self.geo_reference.get_zone() is not DEFAULT_ZONE, msg
383            lats_longs = []
384            for point in self.get_data_points(True):
385                ### UTMtoLL(northing, easting, zone,
386                lat_calced, long_calced = UTMtoLL(point[1],point[0], 
387                                                  zone, isSouthHemisphere)
388                lats_longs.append((lat_calced, long_calced)) # to hash
389            return lats_longs
390           
391        if absolute is True and geo_reference is None:
392            return self.geo_reference.get_absolute(self.data_points)
393        elif geo_reference is not None:
394            return geo_reference.change_points_geo_ref \
395                               (self.data_points, 
396                                self.geo_reference)
397        else:
398            return self.data_points
399       
400   
401    def get_attributes(self, attribute_name=None):
402        """Return values for one named attribute.
403
404        If attribute_name is None, default_attribute_name is used
405        """
406
407        if attribute_name is None:
408            if self.default_attribute_name is not None:
409                attribute_name = self.default_attribute_name
410            else:
411                attribute_name = self.attributes.keys()[0] 
412                # above line takes the first one from keys
413       
414        if self.verbose is True:
415            print 'Using attribute %s' %attribute_name
416            print 'Available attributes: %s' %(self.attributes.keys())       
417
418        msg = 'Attribute name %s does not exist in data set' %attribute_name
419        assert self.attributes.has_key(attribute_name), msg
420
421        return self.attributes[attribute_name]
422
423    def get_all_attributes(self):
424        """
425        Return values for all attributes.
426        The return value is either None or a dictionary (possibly empty).
427        """
428
429        return self.attributes
430
431    def __add__(self, other):
432        """
433        Returns the addition of 2 geospatical objects,
434        objects are concatencated to the end of each other
435           
436        NOTE: doesn't add if objects contain different
437        attributes 
438       
439        Always return absolute points!
440        This alse means, that if you add None to the object,
441        it will be turned into absolute coordinates
442
443        other can be None in which case nothing is added to self.
444        """
445
446
447        # find objects zone and checks if the same
448        geo_ref1 = self.get_geo_reference()
449        zone1 = geo_ref1.get_zone()
450
451
452        if other is not None:
453
454            geo_ref2 = other.get_geo_reference()
455            zone2 = geo_ref2.get_zone()
456
457            geo_ref1.reconcile_zones(geo_ref2)
458
459       
460            new_points = concatenate((self.get_data_points(absolute=True),
461                                      other.get_data_points(absolute=True)),
462                                      axis = 0)       
463
464       
465     
466            # Concatenate attributes if any
467            if self.attributes is None:
468                if other.attributes is not None:
469                    msg = 'Geospatial data must have the same \n'
470                    msg += 'attributes to allow addition.'
471                    raise Exception, msg
472               
473                new_attributes = None
474            else:   
475                new_attributes = {}
476                for x in self.attributes.keys():
477                    if other.attributes.has_key(x):
478
479                        attrib1 = self.attributes[x]
480                        attrib2 = other.attributes[x]
481                        new_attributes[x] = concatenate((attrib1, attrib2))
482
483                    else:
484                        msg = 'Geospatial data must have the same \n'
485                        msg += 'attributes to allow addition.'
486                        raise Exception, msg
487
488
489        else:
490            #other is None:
491           
492            new_points = self.get_data_points(absolute=True)
493            new_attributes = self.attributes
494
495                   
496
497        # Instantiate new data object and return absolute coordinates
498        new_geo_ref = Geo_reference(geo_ref1.get_zone(), 0.0, 0.0)
499        return Geospatial_data(new_points,
500                               new_attributes,
501                               new_geo_ref)
502
503
504    def __radd__(self, other):
505        """Handle cases like None + Geospatial_data(...)
506        """
507
508        return self + other
509
510   
511   
512    ###
513    #  IMPORT/EXPORT POINTS FILES
514    ###
515
516    def import_points_file(self, file_name, delimiter=None, verbose=False):
517        """ load an .txt, .csv or .xya or .pts file
518        Note: will throw an IOError if it can't load the file.
519        Catch these!
520
521        Post condition: self.attributes dictionary has been set
522        """
523       
524        if access(file_name, F_OK) == 0 :
525            msg = 'File %s does not exist or is not accessible' %file_name
526            raise IOError, msg
527       
528        attributes = {}
529        if file_name[-4:]== ".xya":
530            # Maybe not phase-out, so we can load in geo-ref info
531            #msg = 'Text file format is moving to comma seperated .txt files.'
532            #warn(msg, DeprecationWarning)
533            try:
534                if delimiter == None:
535                    try:
536                        fd = open(file_name)
537                        data_points, attributes, geo_reference =\
538                                     _read_xya_file(fd, ',')
539                    except TitleError:
540                        fd.close()
541                        fd = open(file_name)
542                        data_points, attributes, geo_reference =\
543                                     _read_xya_file(fd, ' ')
544                else:
545                    fd = open(file_name)
546                    data_points, attributes, geo_reference =\
547                                 _read_xya_file(fd, delimiter)
548                fd.close()
549            except (IndexError,ValueError,SyntaxError):
550                fd.close()   
551                msg = 'Could not open file %s ' %file_name
552                msg += 'Check the file location.'
553                raise IOError, msg
554            except IOError, e:
555                fd.close() 
556                # Catch this to add an error message
557                msg = 'Could not open file or incorrect file format %s:%s'\
558                      %(file_name, e)
559                raise IOError, msg
560               
561        elif file_name[-4:]== ".pts":
562            try:
563                data_points, attributes, geo_reference =\
564                             _read_pts_file(file_name, verbose)
565            except IOError, e:   
566                msg = 'Could not open file %s ' %file_name
567                raise IOError, msg 
568       
569        elif file_name[-4:]== ".txt" or file_name[-4:]== ".csv":
570            try:
571                data_points, attributes, geo_reference =\
572                             _read_csv_file(file_name, verbose)
573            except IOError, e:
574                # This should only be if a file is not found
575                msg = 'Could not open file %s. ' %file_name
576                msg += 'Check the file location.'
577                raise IOError, msg
578            except SyntaxError, e:
579                # This should only be if there is a format error
580                msg = 'Could not open file %s. \n' %file_name
581                msg += Error_message['IOError']
582                #print "msg", msg
583                raise SyntaxError, msg
584        else:     
585            msg = 'Extension %s is unknown' %file_name[-4:]
586            raise IOError, msg
587#        print'in import data_points', data_points
588#        print'in import attributes', attributes
589#        print'in import data_points', geo_reference
590        self.data_points = data_points
591        self.attributes = attributes
592        self.geo_reference = geo_reference
593   
594#        return all_data
595   
596    def export_points_file(self, file_name, absolute=True, 
597                           as_lat_long=False, isSouthHemisphere=True):
598       
599        """
600        write a points file, file_name, as a text (.xya) or binary (.pts) file
601        file_name is the file name, including the extension
602        The point_dict is defined at the top of this file.
603       
604        If absolute is True data the xll and yll are added to the points value
605        and the xll and yll of the geo_reference are set to 0.
606       
607        If absolute is False data points at returned as relative to the xll
608        and yll and geo_reference remains uneffected
609       
610        isSouthHemisphere: is only used when getting data
611        points "as_lat_long" is True and if FALSE will return lats and
612        longs valid for the Northern Hemisphere.
613       
614        """
615
616        if absolute is False and file_name[-4:] == ".xya":
617            msg = 'The text file values must be absolute.   '
618            msg += 'Text file format is moving to comma seperated .txt files.'
619            warn(msg, DeprecationWarning) 
620
621        if (file_name[-4:] == ".xya"):
622            msg = '.xya format is deprecated.  Please use .txt.'
623            warn(msg, DeprecationWarning)
624            if absolute is True:     
625                geo_ref = deepcopy(self.geo_reference)
626                geo_ref.xllcorner = 0
627                geo_ref.yllcorner = 0   
628                _write_xya_file(file_name,
629                                self.get_data_points(absolute=True), 
630                                self.get_all_attributes(),
631                                geo_ref)
632            else:
633                _write_xya_file(file_name,
634                                self.get_data_points(absolute=False), 
635                                self.get_all_attributes(),
636                                self.get_geo_reference())
637                                   
638        elif (file_name[-4:] == ".pts"):
639            if absolute is True:
640                geo_ref = deepcopy(self.geo_reference)
641                geo_ref.xllcorner = 0
642                geo_ref.yllcorner = 0
643                _write_pts_file(file_name,
644                                self.get_data_points(absolute), 
645                                self.get_all_attributes(),
646                                geo_ref)
647            else:
648                _write_pts_file(file_name,
649                                self.get_data_points(absolute), 
650                                self.get_all_attributes(),
651                                self.get_geo_reference())
652               
653        elif file_name[-4:] == ".txt" or file_name[-4:] == ".csv":
654            msg = "ERROR: trying to write a .txt file with relative data."
655            assert absolute, msg
656            _write_csv_file(file_name,
657                            self.get_data_points(absolute=True,
658                                                 as_lat_long=as_lat_long,
659                                  isSouthHemisphere=isSouthHemisphere), 
660                            self.get_all_attributes(),
661                            as_lat_long=as_lat_long)
662             
663        elif file_name[-4:] == ".urs" :
664            msg = "ERROR: Can not write a .urs file as a relative file."
665            assert absolute, msg
666            _write_urs_file(file_name,
667                            self.get_data_points(as_lat_long=True,
668                                  isSouthHemisphere=isSouthHemisphere))
669                                                         
670        else:
671            msg = 'Unknown file type %s ' %file_name
672            raise IOError, msg
673       
674    def get_sample(self, indices):
675        """ Returns a object which is a subset of the original
676        and the data points and attributes in this new object refer to
677        the indices provided
678       
679        Input
680            indices- a list of integers that represent the new object
681        Output
682            New geospatial data object representing points specified by
683            the indices
684            """
685        #FIXME: add the geo_reference to this
686       
687        points = self.get_data_points()
688        sampled_points = take(points, indices)
689
690        attributes = self.get_all_attributes()
691
692        sampled_attributes = {}
693        if attributes is not None:
694            for key, att in attributes.items():
695                sampled_attributes[key] = take(att, indices)
696
697        return Geospatial_data(sampled_points, sampled_attributes)
698   
699   
700    def split(self, factor=0.5, verbose=False):
701        """Returns two geospatial_data object, first is the size of the 'factor'
702        smaller the original and the second is the remainder. The two new
703        object are disjoin set of each other.
704       
705        Points of the two new object have selected RANDOMLY.
706        AND if factor is a decimal it will round (2.25 to 2 and 2.5 to 3)
707       
708        This method create two lists of indices which are passed into get_sample.
709        The lists are created using random numbers, and they are unique sets eg.
710        total_list(1,2,3,4,5,6,7,8,9)  random_list(1,3,6,7,9) and remainder_list(0,2,4,5,8)
711               
712        Input - the factor which to split the object, if 0.1 then 10% of the
713            object will be returned
714       
715        Output - two geospatial_data objects that are disjoint sets of the
716            original
717            """
718       
719        i=0
720        self_size = len(self)
721        random_list = []
722        remainder_list = []
723        new_size = round(factor*self_size)
724       
725        #find unique random numbers
726        if verbose: print "make unique random number list and get indices"
727
728        total=array(range(self_size))
729        total_list = total.tolist()
730        if verbose: print "total list len",len(total_list)
731               
732        #there will be repeated random numbers however will not be a
733        #problem as they are being 'pop'ed out of array so if there
734        #are two numbers the same they will pop different indicies,
735        #still basically random
736        ## create list of non-unquie random numbers
737        if verbose: print "create random numbers list %s long" %new_size
738        random_num = randint(0,self_size-1,(int(new_size),))
739        random_num = random_num.tolist()
740
741        #need to sort and reverse so the pop() works correctly
742        random_num.sort()
743        random_num.reverse()       
744       
745        if verbose: print "make random number list and get indices"
746        j=0
747        k=1
748        remainder_list = total_list[:]
749        #pops array index (random_num) from remainder_list (which starts as the
750        #total_list and appends to random_list 
751        random_num_len = len(random_num)
752        for i in random_num:
753            random_list.append(remainder_list.pop(i))
754            j+=1
755            #prints progress
756            if verbose and round(random_num_len/10*k)==j:
757                print '(%s/%s)' %(j, random_num_len)
758                k+=1
759       
760        #FIXME: move to tests, it might take a long time
761        #then create an array of random lenght between 500 and 1000,
762        #and use a random factor between 0 and 1
763        #setup for assertion
764        test_total = random_list[:]
765        test_total.extend(remainder_list)
766        test_total.sort() 
767        msg = 'The two random lists made from the original list when added together '\
768         'DO NOT equal the original list'
769        assert (total_list==test_total),msg
770
771        #get new samples
772        if verbose: print "get values of indices for random list"
773        G1 = self.get_sample(random_list)
774        if verbose: print "get values of indices for opposite of random list"
775        G2 = self.get_sample(remainder_list)
776
777        return G1, G2
778
779    def __iter__(self):
780        """Read in the header, number_of_points and save the
781        file pointer position
782        """
783
784        from Scientific.IO.NetCDF import NetCDFFile
785       
786        #FIXME - what to do if the file isn't there
787
788        # FIXME (Ole): Shouldn't this go into the constructor?
789        # ... and shouldn't it be called block_size?
790        if self.max_read_lines is None:
791            self.max_read_lines = MAX_READ_LINES
792       
793
794        if self.file_name[-4:] == ".xya":
795            # FIXME (Ole): shouldn't the xya format be replaced by txt/csv?
796            #  Currently both file formats are used.
797
798            # FIXME (Ole): This has to go - it caused Ted Rigby to waste
799            # time trying to read in his data in xya format with an
800            # inevitable memory error appearing.
801           
802            #let's just read it all
803            msg = 'The xya format is deprecated. Use csv or pts.'
804            warn(msg, DeprecationWarning)             
805           
806        elif self.file_name[-4:] == ".pts":
807           
808            # see if the file is there.  Throw a QUIET IO error if it isn't
809            fd = open(self.file_name,'r')
810            fd.close()
811   
812            #throws prints to screen if file not present
813            self.fid = NetCDFFile(self.file_name, 'r')
814           
815            self.blocking_georef, self.blocking_keys, self.number_of_points =\
816                                  _read_pts_file_header(self.fid,
817                                                        self.verbose)
818            self.start_row = 0
819            self.last_row = self.number_of_points           
820            self.show_verbose = 0
821            self.verbose_block_size = (self.last_row + 10)/10
822            self.block_number = 0
823            self.number_of_blocks = self.number_of_points/self.max_read_lines
824            # This computes the number of full blocks. The last block may be
825            # smaller and won't be ircluded in this estimate.
826           
827            if self.verbose is True:
828                print 'Reading %d points (in ~%d blocks) from file %s. '\
829                      %(self.number_of_points,
830                        self.number_of_blocks,
831                        self.file_name),
832                print 'Each block consists of %d data points'\
833                      %self.max_read_lines
834           
835        else:
836            # assume the file is a csv file
837            file_pointer = open(self.file_name)
838            self.header, self.file_pointer = \
839                         _read_csv_file_header(file_pointer)
840            self.blocking_georef = None # Used for reconciling zones
841
842        return self
843   
844   
845    def next(self):
846        """ read a block, instanciate a new geospatial and return it"""
847       
848        if self.file_name[-4:]== ".xya" :
849            # FIXME (Ole): shouldn't the xya format be replaced by txt/csv?
850            #  Currently both file formats are used.
851           
852            if not hasattr(self,'finished_reading') or \
853                   self.finished_reading is False:
854                #let's just read it all
855                geo = Geospatial_data(self.file_name)
856                self.finished_reading = True
857            else:
858                raise StopIteration
859                self.finished_reading = False
860               
861        elif self.file_name[-4:] == ".pts":
862            if self.start_row == self.last_row:
863                # read the end of the file last iteration
864                # remove blocking attributes
865                self.fid.close()
866                del self.max_read_lines
867                del self.blocking_georef
868                del self.last_row
869                del self.start_row
870                del self.blocking_keys
871                del self.fid
872                raise StopIteration
873            fin_row = self.start_row + self.max_read_lines
874            if fin_row > self.last_row:
875                fin_row = self.last_row
876
877               
878           
879            if self.verbose is True:
880                if self.show_verbose >= self.start_row and \
881                       self.show_verbose < fin_row:
882
883                   
884                    #print 'Doing %d of %d' %(self.start_row, self.last_row+10)
885
886                    print 'Reading block %d (points %d to %d) out of %d'\
887                          %(self.block_number,
888                            self.start_row,
889                            fin_row,
890                            self.number_of_blocks)
891
892                   
893                    self.show_verbose += max(self.max_read_lines,
894                                             self.verbose_block_size)
895
896                 
897            # Read next block
898            pointlist, att_dict, = \
899                       _read_pts_file_blocking(self.fid,
900                                               self.start_row,
901                                               fin_row,
902                                               self.blocking_keys)
903           
904            geo = Geospatial_data(pointlist, att_dict, self.blocking_georef)
905            self.start_row = fin_row
906           
907            self.block_number += 1           
908           
909        else:
910            # assume the file is a csv file
911            try:
912                pointlist, att_dict, geo_ref, self.file_pointer = \
913                   _read_csv_file_blocking( self.file_pointer,
914                                            self.header[:],
915                                            max_read_lines=self.max_read_lines,
916                                            verbose=self.verbose)
917
918                # Check that the zones haven't changed.
919                if geo_ref is not None:
920                    geo_ref.reconcile_zones(self.blocking_georef)
921                    self.blocking_georef = geo_ref
922                elif self.blocking_georef is not None:
923                   
924                    msg = 'Geo reference given, then not given.'
925                    msg += ' This should not happen.' 
926                    raise ValueError, msg
927                geo = Geospatial_data(pointlist, att_dict, geo_ref)
928            except StopIteration:
929                self.file_pointer.close()
930                del self.header
931                del self.file_pointer
932                raise StopIteration
933            except ANUGAError:
934                self.file_pointer.close()
935                del self.header
936                del self.file_pointer
937                raise
938            except SyntaxError:
939                self.file_pointer.close()
940                del self.header
941                del self.file_pointer
942                # This should only be if there is a format error
943                msg = 'Could not open file %s. \n' %self.file_name
944                msg += Error_message['IOError']
945                raise SyntaxError, msg
946        return geo
947##################### Error messages ###########
948Error_message = {}
949Em = Error_message
950Em['IOError'] = "NOTE: The format for a comma separated .txt/.csv file is:\n"
951Em['IOError'] += "        1st line:     [column names]\n"
952Em['IOError'] += "        other lines:  [x value], [y value], [attributes]\n"
953Em['IOError'] += "\n"
954Em['IOError'] += "           for example:\n"
955Em['IOError'] += "           x, y, elevation, friction\n"
956Em['IOError'] += "           0.6, 0.7, 4.9, 0.3\n"
957Em['IOError'] += "           1.9, 2.8, 5, 0.3\n"
958Em['IOError'] += "           2.7, 2.4, 5.2, 0.3\n"
959Em['IOError'] += "\n"
960Em['IOError'] += "The first two columns are assumed to be x, y coordinates.\n"
961Em['IOError'] += "The attribute values must be numeric.\n"
962
963def _set_using_lat_long(latitudes,
964                        longitudes,
965                        geo_reference,
966                        data_points,
967                        points_are_lats_longs):
968    """
969    if the points has lat long info, assume it is in (lat, long) order.
970    """
971   
972    if geo_reference is not None:
973        msg = """A georeference is specified yet latitude and longitude
974        are also specified!"""
975        raise ValueError, msg
976   
977    if data_points is not None and not points_are_lats_longs:
978        msg = """Data points are specified yet latitude and
979        longitude are also specified!"""
980        raise ValueError, msg
981   
982    if points_are_lats_longs:
983        if data_points is None:
984            msg = """Data points are not specified !"""
985            raise ValueError, msg
986        lats_longs = ensure_numeric(data_points)
987        latitudes = ravel(lats_longs[:,0:1])
988        longitudes = ravel(lats_longs[:,1:])
989       
990    if latitudes is None and longitudes is None:
991        msg = """Latitudes and Longitudes are not."""
992        raise ValueError, msg
993   
994    if latitudes is None:
995        msg = """Longitudes are specified yet latitudes aren't!"""
996        raise ValueError, msg
997   
998    if longitudes is None:
999        msg = """Latitudes are specified yet longitudes aren't!"""
1000        raise ValueError, msg
1001   
1002    data_points, zone  = convert_from_latlon_to_utm(latitudes=latitudes,
1003                                                    longitudes=longitudes)
1004    return data_points, Geo_reference(zone=zone)
1005
1006   
1007def _read_pts_file(file_name, verbose=False):
1008    """Read .pts NetCDF file
1009   
1010    Return a dic of array of points, and dic of array of attribute
1011    eg
1012    dic['points'] = [[1.0,2.0],[3.0,5.0]]
1013    dic['attributelist']['elevation'] = [[7.0,5.0]
1014    """   
1015
1016    from Scientific.IO.NetCDF import NetCDFFile
1017   
1018    if verbose: print 'Reading ', file_name
1019   
1020       
1021    # see if the file is there.  Throw a QUIET IO error if it isn't
1022    fd = open(file_name,'r')
1023    fd.close()
1024   
1025    #throws prints to screen if file not present
1026    fid = NetCDFFile(file_name, 'r') 
1027   
1028    pointlist = array(fid.variables['points'])
1029    keys = fid.variables.keys()
1030    if verbose: print 'Got %d variables: %s' %(len(keys), keys)
1031    try:
1032        keys.remove('points')
1033    except IOError, e:       
1034        fid.close()   
1035        msg = 'Expected keyword "points" but could not find it'
1036        raise IOError, msg
1037   
1038    attributes = {}
1039    for key in keys:
1040        if verbose: print "reading attribute '%s'" %key
1041           
1042        attributes[key] = array(fid.variables[key])
1043   
1044   
1045    try:
1046        geo_reference = Geo_reference(NetCDFObject=fid)
1047    except AttributeError, e:
1048        geo_reference = None
1049   
1050    fid.close()
1051   
1052    return pointlist, attributes, geo_reference
1053
1054
1055def _read_csv_file(file_name, verbose=False):
1056    """Read .csv file
1057   
1058    Return a dic of array of points, and dic of array of attribute
1059    eg
1060    dic['points'] = [[1.0,2.0],[3.0,5.0]]
1061    dic['attributelist']['elevation'] = [[7.0,5.0]
1062    """
1063   
1064    file_pointer = open(file_name)
1065    header, file_pointer = _read_csv_file_header(file_pointer)
1066    try:
1067        pointlist, att_dict, geo_ref, file_pointer = \
1068                   _read_csv_file_blocking( \
1069                file_pointer,
1070                header,
1071                max_read_lines=1e30) #If the file is bigger that this, block..
1072    except ANUGAError:
1073        file_pointer.close()
1074        raise 
1075    file_pointer.close()
1076    return pointlist, att_dict, geo_ref   
1077
1078CSV_DELIMITER = ','
1079def _read_csv_file_header(file_pointer,
1080                          delimiter=CSV_DELIMITER,
1081                          verbose=False):
1082
1083    """Read the header of a .csv file
1084    Return a list of the header names
1085    """
1086    line = file_pointer.readline()
1087    header = clean_line(line, delimiter)
1088    return header, file_pointer
1089
1090def _read_csv_file_blocking(file_pointer, header,
1091                            delimiter=CSV_DELIMITER,
1092                            max_read_lines=MAX_READ_LINES,
1093                            verbose=False):
1094   
1095
1096    """
1097    Read the body of a .csv file.
1098    header: The list header of the csv file, with the x and y labels.
1099    """
1100    points = []
1101    pointattributes = []
1102    att_dict = {}
1103
1104    #This is to remove the x and y headers.
1105    header = header[:]
1106    try:
1107        x_header = header.pop(0)
1108        y_header = header.pop(0)
1109    except IndexError:
1110        # if there are not two columns this will occur.
1111        # eg if it is a space seperated file
1112        raise SyntaxError
1113   
1114    read_lines = 0
1115    while read_lines<max_read_lines:
1116        line = file_pointer.readline()
1117        #print "line",line
1118        numbers = clean_line(line,delimiter)
1119        if len(numbers) <= 1:
1120            break
1121        if line[0] == '#':
1122            continue
1123        read_lines += 1
1124        try:
1125            x = float(numbers[0])
1126            y = float(numbers[1])
1127            points.append([x,y])
1128            numbers.pop(0)
1129            numbers.pop(0)
1130            if len(header) != len(numbers):
1131                file_pointer.close() 
1132                msg = "File load error.  There might be a problem with the file header"
1133                raise SyntaxError, msg
1134            for i,num in enumerate(numbers):
1135                num.strip()
1136                if num != '\n' and num != '':
1137                    #attributes.append(float(num))
1138                    att_dict.setdefault(header[i],[]).append(float(num))
1139        #except IOError:           
1140        except ValueError:
1141            raise SyntaxError
1142    if points == []:
1143        raise StopIteration
1144       
1145   
1146    pointlist = array(points).astype(Float)
1147    for key in att_dict.keys():
1148        att_dict[key] = array(att_dict[key]).astype(Float)
1149
1150    #Do stuff here so the info is in lat's and longs
1151    geo_ref = None
1152    x_header = lower(x_header[:3])
1153    y_header = lower(y_header[:3])
1154    if (x_header == 'lon' or  x_header == 'lat') and \
1155       (y_header == 'lon' or  y_header == 'lat'):
1156        if x_header == 'lon':
1157            longitudes = ravel(pointlist[:,0:1])
1158            latitudes = ravel(pointlist[:,1:])
1159        else:
1160            latitudes = ravel(pointlist[:,0:1])
1161            longitudes = ravel(pointlist[:,1:])
1162       
1163        pointlist, geo_ref = _set_using_lat_long(latitudes,
1164                                                 longitudes,
1165                                                 geo_reference=None,
1166                                                 data_points=None,
1167                                                 points_are_lats_longs=False)
1168    return pointlist, att_dict, geo_ref, file_pointer
1169
1170def _read_pts_file_header(fid, verbose=False):
1171
1172    """Read the geo_reference and number_of_points from a .pts file
1173    """
1174   
1175    keys = fid.variables.keys()
1176    try:
1177        keys.remove('points')
1178    except IOError, e:       
1179        fid.close()   
1180        msg = 'Expected keyword "points" but could not find it'
1181        raise IOError, msg
1182    if verbose: print 'Got %d variables: %s' %(len(keys), keys)
1183   
1184    try:
1185        geo_reference = Geo_reference(NetCDFObject=fid)
1186    except AttributeError, e:
1187        geo_reference = None
1188
1189    return geo_reference, keys, fid.dimensions['number_of_points']
1190
1191def _read_pts_file_blocking(fid, start_row, fin_row, keys):
1192                            #verbose=False):
1193   
1194
1195    """
1196    Read the body of a .csv file.
1197    header: The list header of the csv file, with the x and y labels.
1198    """
1199   
1200    pointlist = array(fid.variables['points'][start_row:fin_row])
1201   
1202    attributes = {}
1203    for key in keys:
1204        attributes[key] = array(fid.variables[key][start_row:fin_row])
1205
1206    return pointlist, attributes
1207   
1208   
1209def _read_xya_file(fd, delimiter):
1210    points = []
1211    pointattributes = []
1212    title = fd.readline()
1213    att_names = clean_line(title,delimiter)
1214    att_dict = {}
1215    line = fd.readline()
1216    numbers = clean_line(line,delimiter)
1217   
1218    while len(numbers) > 1 and line[0] <> '#':
1219        if numbers != []:
1220            try:
1221                x = float(numbers[0])
1222                y = float(numbers[1])
1223                points.append([x,y])
1224                numbers.pop(0)
1225                numbers.pop(0)
1226                if len(att_names) != len(numbers):
1227                    fd.close()
1228                    # It might not be a problem with the title
1229                    #raise TitleAmountError
1230                    raise IOError
1231                for i,num in enumerate(numbers):
1232                    num.strip()
1233                    if num != '\n' and num != '':
1234                        #attributes.append(float(num))
1235                        att_dict.setdefault(att_names[i],[]).append(float(num))
1236            except ValueError:
1237                raise SyntaxError
1238        line = fd.readline()
1239        numbers = clean_line(line,delimiter)
1240   
1241    if line == '':
1242        geo_reference = None
1243    else:
1244        geo_reference = Geo_reference(ASCIIFile=fd,read_title=line)
1245       
1246   
1247    pointlist = array(points).astype(Float)
1248    for key in att_dict.keys():
1249        att_dict[key] = array(att_dict[key]).astype(Float)
1250   
1251    return pointlist, att_dict, geo_reference
1252
1253def _write_pts_file(file_name,
1254                    write_data_points,
1255                    write_attributes=None, 
1256                    write_geo_reference=None):
1257    """
1258    Write .pts NetCDF file   
1259
1260    NOTE: Below might not be valid ask Duncan : NB 5/2006
1261   
1262    WARNING: This function mangles the point_atts data structure
1263    #F??ME: (DSG)This format has issues.
1264    # There can't be an attribute called points
1265    # consider format change
1266    # method changed by NB not sure if above statement is correct
1267   
1268    should create new test for this
1269    legal_keys = ['pointlist', 'attributelist', 'geo_reference']
1270    for key in point_atts.keys():
1271        msg = 'Key %s is illegal. Valid keys are %s' %(key, legal_keys)
1272        assert key in legal_keys, msg
1273    """   
1274    from Scientific.IO.NetCDF import NetCDFFile
1275    # NetCDF file definition
1276    outfile = NetCDFFile(file_name, 'w')
1277   
1278    #Create new file
1279    outfile.institution = 'Geoscience Australia'
1280    outfile.description = 'NetCDF format for compact and portable storage ' +\
1281                          'of spatial point data'
1282   
1283    # dimension definitions
1284    shape = write_data_points.shape[0]
1285    outfile.createDimension('number_of_points', shape) 
1286    outfile.createDimension('number_of_dimensions', 2) #This is 2d data
1287   
1288    # variable definition
1289    outfile.createVariable('points', Float, ('number_of_points',
1290                                             'number_of_dimensions'))
1291
1292    #create variables 
1293    outfile.variables['points'][:] = write_data_points #.astype(Float32)
1294
1295    if write_attributes is not None:
1296        for key in write_attributes.keys():
1297            outfile.createVariable(key, Float, ('number_of_points',))
1298            outfile.variables[key][:] = write_attributes[key] #.astype(Float32)
1299       
1300    if write_geo_reference is not None:
1301        write_NetCDF_georeference(write_geo_reference, outfile)
1302       
1303    outfile.close() 
1304 
1305
1306
1307def _write_xya_file(file_name,
1308                    write_data_points,
1309                    write_attributes=None, 
1310                    write_geo_reference=None, 
1311                    delimiter=','):
1312    """
1313    export a file, file_name, with the xya format
1314   
1315    """
1316    points = write_data_points
1317    pointattributes = write_attributes
1318   
1319    fd = open(file_name,'w')
1320    titlelist = ""
1321    if pointattributes is not None:   
1322        for title in pointattributes.keys():
1323            titlelist = titlelist + title + delimiter
1324        titlelist = titlelist[0:-len(delimiter)] # remove the last delimiter
1325    fd.write(titlelist+"\n")
1326   
1327    #<vertex #> <x> <y> [attributes]
1328    for i, vert in enumerate( points):
1329
1330
1331        if pointattributes is not None:           
1332            attlist = ","
1333            for att in pointattributes.keys():
1334                attlist = attlist + str(pointattributes[att][i])+ delimiter
1335            attlist = attlist[0:-len(delimiter)] # remove the last delimiter
1336            attlist.strip()
1337        else:
1338            attlist = ''
1339
1340        fd.write(str(vert[0]) + delimiter +
1341                 str(vert[1]) + attlist + "\n")
1342
1343    if  write_geo_reference is not None:
1344        write_geo_reference = ensure_geo_reference(write_geo_reference)
1345        write_geo_reference.write_ASCII(fd)
1346    fd.close()
1347
1348
1349def _write_csv_file(file_name,
1350                    write_data_points,
1351                    write_attributes=None,
1352                    as_lat_long=False,
1353                    delimiter=','):
1354    """
1355    export a file, file_name, with the csv format
1356   
1357    """
1358    points = write_data_points
1359    pointattributes = write_attributes
1360   
1361    fd = open(file_name,'w')
1362    if as_lat_long:
1363        titlelist = "latitude" + delimiter + "longitude"  + delimiter
1364    else:
1365        titlelist = "x" + delimiter + "y"  + delimiter
1366    if pointattributes is not None:   
1367        for title in pointattributes.keys():
1368            titlelist = titlelist + title + delimiter
1369        titlelist = titlelist[0:-len(delimiter)] # remove the last delimiter
1370    fd.write(titlelist+"\n")
1371   
1372    # <x/lat> <y/long> [attributes]
1373    for i, vert in enumerate( points):
1374
1375
1376        if pointattributes is not None:           
1377            attlist = ","
1378            for att in pointattributes.keys():
1379                attlist = attlist + str(pointattributes[att][i])+ delimiter
1380            attlist = attlist[0:-len(delimiter)] # remove the last delimiter
1381            attlist.strip()
1382        else:
1383            attlist = ''
1384
1385        fd.write(str(vert[0]) + delimiter +
1386                 str(vert[1]) + attlist + "\n")
1387
1388    fd.close()
1389 
1390def _write_urs_file(file_name,
1391                    points,
1392                    delimiter=' '):
1393    """
1394    export a file, file_name, with the urs format
1395    the data points are in lats and longs
1396   
1397    """   
1398    fd = open(file_name,'w')   
1399    fd.write(str(len(points))+"\n")   
1400    # <lat> <long> <id#>
1401    for i, vert in enumerate( points):
1402        fd.write(str(round(vert[0],7)) + delimiter + \
1403                 str(round(vert[1],7)) + delimiter +str(i)+ "\n")
1404    fd.close()
1405   
1406def _point_atts2array(point_atts):
1407    point_atts['pointlist'] = array(point_atts['pointlist']).astype(Float)
1408   
1409    for key in point_atts['attributelist'].keys():
1410        point_atts['attributelist'][key]=\
1411                array(point_atts['attributelist'][key]).astype(Float)
1412    return point_atts
1413
1414   
1415
1416
1417def geospatial_data2points_dictionary(geospatial_data):
1418    """Convert geospatial data to points_dictionary
1419    """
1420
1421    points_dictionary = {}
1422    points_dictionary['pointlist'] = geospatial_data.data_points
1423
1424    points_dictionary['attributelist'] = {}
1425
1426    for attribute_name in geospatial_data.attributes.keys():
1427        val = geospatial_data.attributes[attribute_name]
1428        points_dictionary['attributelist'][attribute_name] = val
1429
1430    points_dictionary['geo_reference'] = geospatial_data.geo_reference
1431
1432    return points_dictionary
1433
1434   
1435def points_dictionary2geospatial_data(points_dictionary):
1436    """Convert points_dictionary to geospatial data object
1437    """
1438
1439    msg = 'Points dictionary must have key pointlist' 
1440    assert points_dictionary.has_key('pointlist'), msg
1441
1442    msg = 'Points dictionary must have key attributelist'     
1443    assert points_dictionary.has_key('attributelist'), msg       
1444
1445    if points_dictionary.has_key('geo_reference'):
1446        geo = points_dictionary['geo_reference']
1447    else:
1448        geo = None
1449   
1450    return Geospatial_data(points_dictionary['pointlist'],
1451                           points_dictionary['attributelist'],
1452                           geo_reference = geo)
1453
1454def clean_line(line,delimiter):     
1455    """Remove whitespace
1456    """
1457    #print ">%s" %line
1458    line = line.strip()
1459    #print "stripped>%s" %line
1460    numbers = line.split(delimiter)
1461    i = len(numbers) - 1
1462    while i >= 0:
1463        if numbers[i] == '':
1464            numbers.pop(i)
1465        else:
1466            numbers[i] = numbers[i].strip()
1467       
1468        i += -1
1469    #for num in numbers:
1470    #    print "num>%s<" %num
1471    return numbers
1472           
1473def ensure_absolute(points, geo_reference=None):
1474    """
1475    This function inputs several formats and
1476    outputs one format. - a numeric array of absolute points.
1477
1478    Inputed formats are;
1479    points: List or numeric array of coordinate pairs [xi, eta] of
1480              points or geospatial object or points file name
1481
1482    mesh_origin: A geo_reference object or 3-tuples consisting of
1483                 UTM zone, easting and northing.
1484                 If specified vertex coordinates are assumed to be
1485                 relative to their respective origins.
1486    """
1487    if isinstance(points,type('')):
1488        #It's a string
1489        #assume it is a point file
1490        points = Geospatial_data(file_name = points)
1491       
1492    if isinstance(points,Geospatial_data):
1493        points = points.get_data_points( \
1494                absolute = True)
1495        msg = "Use a Geospatial_data object or a mesh origin. Not both."
1496        assert geo_reference == None, msg
1497           
1498    else:
1499        points = ensure_numeric(points, Float)
1500    if geo_reference is None:
1501        geo = None #Geo_reference()
1502    else:
1503        if isinstance(geo_reference, Geo_reference):
1504            geo = geo_reference
1505        else:
1506            geo = Geo_reference(geo_reference[0],
1507                                geo_reference[1],
1508                                geo_reference[2])
1509        points = geo.get_absolute(points)
1510    return points
1511     
1512
1513def ensure_geospatial(points, geo_reference=None):
1514    """
1515    This function inputs several formats and
1516    outputs one format. - a geospatial_data instance.
1517
1518    Inputed formats are;
1519    points: List or numeric array of coordinate pairs [xi, eta] of
1520              points or geospatial object
1521
1522    mesh_origin: A geo_reference object or 3-tuples consisting of
1523                 UTM zone, easting and northing.
1524                 If specified vertex coordinates are assumed to be
1525                 relative to their respective origins.
1526    """
1527    if isinstance(points,Geospatial_data):
1528        msg = "Use a Geospatial_data object or a mesh origin. Not both."
1529        assert geo_reference == None, msg
1530        return points   
1531    else:
1532        points = ensure_numeric(points, Float)
1533    if geo_reference is None:
1534        geo = None
1535    else:
1536        if isinstance(geo_reference, Geo_reference):
1537            geo = geo_reference
1538        else:
1539            geo = Geo_reference(geo_reference[0],
1540                                geo_reference[1],
1541                                geo_reference[2])
1542    points = Geospatial_data(data_points=points, geo_reference=geo)       
1543    return points
1544
1545#def file2xya(filename):
1546   
1547#    G = Geospatial_data(filename)
1548#    G.export_points_file()
1549
1550             
1551
1552         
1553if __name__ == "__main__":
1554    pass
1555   
Note: See TracBrowser for help on using the repository browser.