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

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

Improved verbose reporting when large data files are
processed block wise. This works best for pts files,
but would be good to do for txt and cvs as well.

File size: 54.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):
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        Default: absolute is True.
373        """
374        if as_lat_long is True:
375            msg = "Points need a zone to be converted into lats and longs"
376            assert self.geo_reference is not None, msg
377            zone = self.geo_reference.get_zone()
378            assert self.geo_reference.get_zone() is not DEFAULT_ZONE, msg
379            lats_longs = []
380            for point in self.get_data_points(True):
381                ### UTMtoLL(northing, easting, zone,
382                lat_calced, long_calced = UTMtoLL(point[1],point[0], zone)
383                lats_longs.append((lat_calced, long_calced)) # to hash
384            return lats_longs
385           
386        if absolute is True and geo_reference is None:
387            return self.geo_reference.get_absolute(self.data_points)
388        elif geo_reference is not None:
389            return geo_reference.change_points_geo_ref \
390                               (self.data_points, 
391                                self.geo_reference)
392        else:
393            return self.data_points
394       
395   
396    def get_attributes(self, attribute_name=None):
397        """Return values for one named attribute.
398
399        If attribute_name is None, default_attribute_name is used
400        """
401
402        if attribute_name is None:
403            if self.default_attribute_name is not None:
404                attribute_name = self.default_attribute_name
405            else:
406                attribute_name = self.attributes.keys()[0] 
407                # above line takes the first one from keys
408       
409        if self.verbose is True:
410            print 'Using attribute %s' %attribute_name
411            print 'Available attributes: %s' %(self.attributes.keys())       
412
413        msg = 'Attribute name %s does not exist in data set' %attribute_name
414        assert self.attributes.has_key(attribute_name), msg
415
416        return self.attributes[attribute_name]
417
418    def get_all_attributes(self):
419        """
420        Return values for all attributes.
421        The return value is either None or a dictionary (possibly empty).
422        """
423
424        return self.attributes
425
426    def __add__(self, other):
427        """
428        Returns the addition of 2 geospatical objects,
429        objects are concatencated to the end of each other
430           
431        NOTE: doesn't add if objects contain different
432        attributes 
433       
434        Always return absolute points!
435        This alse means, that if you add None to the object,
436        it will be turned into absolute coordinates
437
438        other can be None in which case nothing is added to self.
439        """
440
441
442        # find objects zone and checks if the same
443        geo_ref1 = self.get_geo_reference()
444        zone1 = geo_ref1.get_zone()
445
446
447        if other is not None:
448
449            geo_ref2 = other.get_geo_reference()
450            zone2 = geo_ref2.get_zone()
451
452            geo_ref1.reconcile_zones(geo_ref2)
453
454       
455            new_points = concatenate((self.get_data_points(absolute=True),
456                                      other.get_data_points(absolute=True)),
457                                      axis = 0)       
458
459       
460     
461            # Concatenate attributes if any
462            if self.attributes is None:
463                if other.attributes is not None:
464                    msg = 'Geospatial data must have the same \n'
465                    msg += 'attributes to allow addition.'
466                    raise Exception, msg
467               
468                new_attributes = None
469            else:   
470                new_attributes = {}
471                for x in self.attributes.keys():
472                    if other.attributes.has_key(x):
473
474                        attrib1 = self.attributes[x]
475                        attrib2 = other.attributes[x]
476                        new_attributes[x] = concatenate((attrib1, attrib2))
477
478                    else:
479                        msg = 'Geospatial data must have the same \n'
480                        msg += 'attributes to allow addition.'
481                        raise Exception, msg
482
483
484        else:
485            #other is None:
486           
487            new_points = self.get_data_points(absolute=True)
488            new_attributes = self.attributes
489
490                   
491
492        # Instantiate new data object and return absolute coordinates
493        new_geo_ref = Geo_reference(geo_ref1.get_zone(), 0.0, 0.0)
494        return Geospatial_data(new_points,
495                               new_attributes,
496                               new_geo_ref)
497
498
499    def __radd__(self, other):
500        """Handle cases like None + Geospatial_data(...)
501        """
502
503        return self + other
504
505   
506   
507    ###
508    #  IMPORT/EXPORT POINTS FILES
509    ###
510
511    def import_points_file(self, file_name, delimiter=None, verbose=False):
512        """ load an .txt, .csv or .xya or .pts file
513        Note: will throw an IOError if it can't load the file.
514        Catch these!
515
516        Post condition: self.attributes dictionary has been set
517        """
518       
519        if access(file_name, F_OK) == 0 :
520            msg = 'File %s does not exist or is not accessible' %file_name
521            raise IOError, msg
522       
523        attributes = {}
524        if file_name[-4:]== ".xya":
525            # Maybe not phase-out, so we can load in geo-ref info
526            #msg = 'Text file format is moving to comma seperated .txt files.'
527            #warn(msg, DeprecationWarning)
528            try:
529                if delimiter == None:
530                    try:
531                        fd = open(file_name)
532                        data_points, attributes, geo_reference =\
533                                     _read_xya_file(fd, ',')
534                    except TitleError:
535                        fd.close()
536                        fd = open(file_name)
537                        data_points, attributes, geo_reference =\
538                                     _read_xya_file(fd, ' ')
539                else:
540                    fd = open(file_name)
541                    data_points, attributes, geo_reference =\
542                                 _read_xya_file(fd, delimiter)
543                fd.close()
544            except (IndexError,ValueError,SyntaxError):
545                fd.close()   
546                msg = 'Could not open file %s ' %file_name
547                msg += 'Check the file location.'
548                raise IOError, msg
549            except IOError, e:
550                fd.close() 
551                # Catch this to add an error message
552                msg = 'Could not open file or incorrect file format %s:%s'\
553                      %(file_name, e)
554                raise IOError, msg
555               
556        elif file_name[-4:]== ".pts":
557            try:
558                data_points, attributes, geo_reference =\
559                             _read_pts_file(file_name, verbose)
560            except IOError, e:   
561                msg = 'Could not open file %s ' %file_name
562                raise IOError, msg 
563       
564        elif file_name[-4:]== ".txt" or file_name[-4:]== ".csv":
565            try:
566                data_points, attributes, geo_reference =\
567                             _read_csv_file(file_name, verbose)
568            except IOError, e:
569                # This should only be if a file is not found
570                msg = 'Could not open file %s. ' %file_name
571                msg += 'Check the file location.'
572                raise IOError, msg
573            except SyntaxError, e:
574                # This should only be if there is a format error
575                msg = 'Could not open file %s. \n' %file_name
576                msg += Error_message['IOError']
577                #print "msg", msg
578                raise SyntaxError, msg
579        else:     
580            msg = 'Extension %s is unknown' %file_name[-4:]
581            raise IOError, msg
582#        print'in import data_points', data_points
583#        print'in import attributes', attributes
584#        print'in import data_points', geo_reference
585        self.data_points = data_points
586        self.attributes = attributes
587        self.geo_reference = geo_reference
588   
589#        return all_data
590   
591    def export_points_file(self, file_name, absolute=True, as_lat_long=False):
592       
593        """
594        write a points file, file_name, as a text (.xya) or binary (.pts) file
595        file_name is the file name, including the extension
596        The point_dict is defined at the top of this file.
597       
598        If absolute is True data the xll and yll are added to the points value
599        and the xll and yll of the geo_reference are set to 0.
600       
601        If absolute is False data points at returned as relative to the xll
602        and yll and geo_reference remains uneffected
603        """
604
605        if absolute is False and file_name[-4:] == ".xya":
606            msg = 'The text file values must be absolute.   '
607            msg += 'Text file format is moving to comma seperated .txt files.'
608            warn(msg, DeprecationWarning) 
609
610        if (file_name[-4:] == ".xya"):
611            msg = '.xya format is deprecated.  Please use .txt.'
612            warn(msg, DeprecationWarning)
613            if absolute is True:     
614                geo_ref = deepcopy(self.geo_reference)
615                geo_ref.xllcorner = 0
616                geo_ref.yllcorner = 0   
617                _write_xya_file(file_name,
618                                self.get_data_points(absolute=True), 
619                                self.get_all_attributes(),
620                                geo_ref)
621            else:
622                _write_xya_file(file_name,
623                                self.get_data_points(absolute=False), 
624                                self.get_all_attributes(),
625                                self.get_geo_reference())
626                                   
627        elif (file_name[-4:] == ".pts"):
628            if absolute is True:
629                geo_ref = deepcopy(self.geo_reference)
630                geo_ref.xllcorner = 0
631                geo_ref.yllcorner = 0
632                _write_pts_file(file_name,
633                                self.get_data_points(absolute), 
634                                self.get_all_attributes(),
635                                geo_ref)
636            else:
637                _write_pts_file(file_name,
638                                self.get_data_points(absolute), 
639                                self.get_all_attributes(),
640                                self.get_geo_reference())
641               
642        elif file_name[-4:] == ".txt" or file_name[-4:] == ".csv":
643            msg = "ERROR: trying to write a .txt file with relative data."
644            assert absolute, msg
645            _write_csv_file(file_name,
646                            self.get_data_points(absolute=True,
647                                                 as_lat_long=as_lat_long), 
648                            self.get_all_attributes(),
649                            as_lat_long=as_lat_long)
650             
651        elif file_name[-4:] == ".urs" :
652            msg = "ERROR: Can not write a .urs file as a relative file."
653            assert absolute, msg
654            _write_urs_file(file_name,
655                            self.get_data_points(as_lat_long=True))
656                                                         
657        else:
658            msg = 'Unknown file type %s ' %file_name
659            raise IOError, msg
660       
661    def get_sample(self, indices):
662        """ Returns a object which is a subset of the original
663        and the data points and attributes in this new object refer to
664        the indices provided
665       
666        Input
667            indices- a list of integers that represent the new object
668        Output
669            New geospatial data object representing points specified by
670            the indices
671            """
672        #FIXME: add the geo_reference to this
673       
674        points = self.get_data_points()
675        sampled_points = take(points, indices)
676
677        attributes = self.get_all_attributes()
678
679        sampled_attributes = {}
680        if attributes is not None:
681            for key, att in attributes.items():
682                sampled_attributes[key] = take(att, indices)
683
684        return Geospatial_data(sampled_points, sampled_attributes)
685   
686   
687    def split(self, factor=0.5, verbose=False):
688        """Returns two geospatial_data object, first is the size of the 'factor'
689        smaller the original and the second is the remainder. The two new
690        object are disjoin set of each other.
691       
692        Points of the two new object have selected RANDOMLY.
693        AND if factor is a decimal it will round (2.25 to 2 and 2.5 to 3)
694       
695        This method create two lists of indices which are passed into get_sample.
696        The lists are created using random numbers, and they are unique sets eg.
697        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)
698               
699        Input - the factor which to split the object, if 0.1 then 10% of the
700            object will be returned
701       
702        Output - two geospatial_data objects that are disjoint sets of the
703            original
704            """
705       
706        i=0
707        self_size = len(self)
708        random_list = []
709        remainder_list = []
710        new_size = round(factor*self_size)
711       
712        #find unique random numbers
713        if verbose: print "make unique random number list and get indices"
714
715        total=array(range(self_size))
716        total_list = total.tolist()
717        if verbose: print "total list len",len(total_list)
718               
719        #there will be repeated random numbers however will not be a
720        #problem as they are being 'pop'ed out of array so if there
721        #are two numbers the same they will pop different indicies,
722        #still basically random
723        ## create list of non-unquie random numbers
724        if verbose: print "create random numbers list %s long" %new_size
725        random_num = randint(0,self_size-1,(int(new_size),))
726        random_num = random_num.tolist()
727
728        #need to sort and reverse so the pop() works correctly
729        random_num.sort()
730        random_num.reverse()       
731       
732        if verbose: print "make random number list and get indices"
733        j=0
734        k=1
735        remainder_list = total_list[:]
736        #pops array index (random_num) from remainder_list (which starts as the
737        #total_list and appends to random_list 
738        random_num_len = len(random_num)
739        for i in random_num:
740            random_list.append(remainder_list.pop(i))
741            j+=1
742            #prints progress
743            if verbose and round(random_num_len/10*k)==j:
744                print '(%s/%s)' %(j, random_num_len)
745                k+=1
746       
747        #FIXME: move to tests, it might take a long time
748        #then create an array of random lenght between 500 and 1000,
749        #and use a random factor between 0 and 1
750        #setup for assertion
751        test_total = random_list[:]
752        test_total.extend(remainder_list)
753        test_total.sort() 
754        msg = 'The two random lists made from the original list when added together '\
755         'DO NOT equal the original list'
756        assert (total_list==test_total),msg
757
758        #get new samples
759        if verbose: print "get values of indices for random list"
760        G1 = self.get_sample(random_list)
761        if verbose: print "get values of indices for opposite of random list"
762        G2 = self.get_sample(remainder_list)
763
764        return G1, G2
765
766    def __iter__(self):
767        """Read in the header, number_of_points and save the
768        file pointer position
769        """
770
771        from Scientific.IO.NetCDF import NetCDFFile
772       
773        #FIXME - what to do if the file isn't there
774
775        # FIXME (Ole): Shouldn't this go into the constructor?
776        # ... and shouldn't it be called block_size?
777        if self.max_read_lines is None:
778            self.max_read_lines = MAX_READ_LINES
779       
780
781        if self.file_name[-4:] == ".xya":
782            # FIXME (Ole): shouldn't the xya format be replaced by txt/csv?
783            #  Currently both file formats are used.
784           
785            #let's just read it all
786            pass
787        elif self.file_name[-4:] == ".pts":
788           
789            # see if the file is there.  Throw a QUIET IO error if it isn't
790            fd = open(self.file_name,'r')
791            fd.close()
792   
793            #throws prints to screen if file not present
794            self.fid = NetCDFFile(self.file_name, 'r')
795           
796            self.blocking_georef, self.blocking_keys, self.number_of_points =\
797                                  _read_pts_file_header(self.fid,
798                                                        self.verbose)
799            self.start_row = 0
800            self.last_row = self.number_of_points           
801            self.show_verbose = 0
802            self.verbose_block_size = (self.last_row + 10)/10
803            self.block_number = 0
804            self.number_of_blocks = self.number_of_points/self.max_read_lines
805            # This computes the number of full blocks. The last block may be
806            # smaller and won't be ircluded in this estimate.
807           
808            if self.verbose is True:
809                print 'Reading %d points (in ~%d blocks) from file %s. '\
810                      %(self.number_of_points,
811                        self.number_of_blocks,
812                        self.file_name),
813                print 'Each block consists of %d data points'\
814                      %self.max_read_lines
815           
816        else:
817            # assume the file is a csv file
818            file_pointer = open(self.file_name)
819            self.header, self.file_pointer = \
820                         _read_csv_file_header(file_pointer)
821            self.blocking_georef = None # Used for reconciling zones
822
823        return self
824   
825   
826    def next(self):
827        """ read a block, instanciate a new geospatial and return it"""
828       
829        if self.file_name[-4:]== ".xya" :
830            # FIXME (Ole): shouldn't the xya format be replaced by txt/csv?
831            #  Currently both file formats are used.
832           
833            if not hasattr(self,'finished_reading') or \
834                   self.finished_reading is False:
835                #let's just read it all
836                geo = Geospatial_data(self.file_name)
837                self.finished_reading = True
838            else:
839                raise StopIteration
840                self.finished_reading = False
841               
842        elif self.file_name[-4:]== ".pts":
843            if self.start_row == self.last_row:
844                # read the end of the file last iteration
845                # remove blocking attributes
846                self.fid.close()
847                del self.max_read_lines
848                del self.blocking_georef
849                del self.last_row
850                del self.start_row
851                del self.blocking_keys
852                del self.fid
853                raise StopIteration
854            fin_row = self.start_row + self.max_read_lines
855            if fin_row > self.last_row:
856                fin_row = self.last_row
857
858               
859           
860            if self.verbose is True:
861                if self.show_verbose >= self.start_row and \
862                       self.show_verbose < fin_row:
863
864                   
865                    #print 'Doing %d of %d' %(self.start_row, self.last_row+10)
866
867                    print 'Reading block %d (points %d to %d) out of %d'\
868                          %(self.block_number,
869                            self.start_row,
870                            fin_row,
871                            self.number_of_blocks)
872
873                   
874                    self.show_verbose += max(self.max_read_lines,
875                                             self.verbose_block_size)
876
877                 
878            # Read next block
879            pointlist, att_dict, = \
880                       _read_pts_file_blocking(self.fid,
881                                               self.start_row,
882                                               fin_row,
883                                               self.blocking_keys)
884           
885            geo = Geospatial_data(pointlist, att_dict, self.blocking_georef)
886            self.start_row = fin_row
887           
888            self.block_number += 1           
889           
890        else:
891            # assume the file is a csv file
892            try:
893                pointlist, att_dict, geo_ref, self.file_pointer = \
894                   _read_csv_file_blocking( self.file_pointer,
895                                            self.header[:],
896                                            max_read_lines=self.max_read_lines,
897                                            verbose=self.verbose)
898
899                # Check that the zones haven't changed.
900                if geo_ref is not None:
901                    geo_ref.reconcile_zones(self.blocking_georef)
902                    self.blocking_georef = geo_ref
903                elif self.blocking_georef is not None:
904                   
905                    msg = 'Geo reference given, then not given.'
906                    msg += ' This should not happen.' 
907                    raise ValueError, msg
908                geo = Geospatial_data(pointlist, att_dict, geo_ref)
909            except StopIteration:
910                self.file_pointer.close()
911                del self.header
912                del self.file_pointer
913                raise StopIteration
914            except ANUGAError:
915                self.file_pointer.close()
916                del self.header
917                del self.file_pointer
918                raise
919            except SyntaxError:
920                self.file_pointer.close()
921                del self.header
922                del self.file_pointer
923                # This should only be if there is a format error
924                msg = 'Could not open file %s. \n' %self.file_name
925                msg += Error_message['IOError']
926                raise SyntaxError, msg
927        return geo
928##################### Error messages ###########
929Error_message = {}
930Em = Error_message
931Em['IOError'] = "NOTE: The format for a comma seperated .txt/.csv file is:\n"
932Em['IOError'] += "        1st line:     [column names]\n"
933Em['IOError'] += "        other lines:  [x value], [y value], [attributes]\n"
934Em['IOError'] += "\n"
935Em['IOError'] += "           for example:\n"
936Em['IOError'] += "           x, y, elevation, friction\n"
937Em['IOError'] += "           0.6, 0.7, 4.9, 0.3\n"
938Em['IOError'] += "           1.9, 2.8, 5, 0.3\n"
939Em['IOError'] += "           2.7, 2.4, 5.2, 0.3\n"
940Em['IOError'] += "\n"
941Em['IOError'] += "The first two columns are assumed to be x, y coordinates.\n"
942Em['IOError'] += "The attribute values must be numeric.\n"
943
944def _set_using_lat_long(latitudes,
945                        longitudes,
946                        geo_reference,
947                        data_points,
948                        points_are_lats_longs):
949    """
950    if the points has lat long info, assume it is in (lat, long) order.
951    """
952   
953    if geo_reference is not None:
954        msg = """A georeference is specified yet latitude and longitude
955        are also specified!"""
956        raise ValueError, msg
957   
958    if data_points is not None and not points_are_lats_longs:
959        msg = """Data points are specified yet latitude and
960        longitude are also specified!"""
961        raise ValueError, msg
962   
963    if points_are_lats_longs:
964        if data_points is None:
965            msg = """Data points are not specified !"""
966            raise ValueError, msg
967        lats_longs = ensure_numeric(data_points)
968        latitudes = ravel(lats_longs[:,0:1])
969        longitudes = ravel(lats_longs[:,1:])
970       
971    if latitudes is None and longitudes is None:
972        msg = """Latitudes and Longitudes are not."""
973        raise ValueError, msg
974   
975    if latitudes is None:
976        msg = """Longitudes are specified yet latitudes aren't!"""
977        raise ValueError, msg
978   
979    if longitudes is None:
980        msg = """Latitudes are specified yet longitudes aren't!"""
981        raise ValueError, msg
982   
983    data_points, zone  = convert_from_latlon_to_utm(latitudes=latitudes,
984                                                    longitudes=longitudes)
985    return data_points, Geo_reference(zone=zone)
986
987   
988def _read_pts_file(file_name, verbose=False):
989    """Read .pts NetCDF file
990   
991    Return a dic of array of points, and dic of array of attribute
992    eg
993    dic['points'] = [[1.0,2.0],[3.0,5.0]]
994    dic['attributelist']['elevation'] = [[7.0,5.0]
995    """   
996
997    from Scientific.IO.NetCDF import NetCDFFile
998   
999    if verbose: print 'Reading ', file_name
1000   
1001       
1002    # see if the file is there.  Throw a QUIET IO error if it isn't
1003    fd = open(file_name,'r')
1004    fd.close()
1005   
1006    #throws prints to screen if file not present
1007    fid = NetCDFFile(file_name, 'r') 
1008   
1009    pointlist = array(fid.variables['points'])
1010    keys = fid.variables.keys()
1011    if verbose: print 'Got %d variables: %s' %(len(keys), keys)
1012    try:
1013        keys.remove('points')
1014    except IOError, e:       
1015        fid.close()   
1016        msg = 'Expected keyword "points" but could not find it'
1017        raise IOError, msg
1018   
1019    attributes = {}
1020    for key in keys:
1021        if verbose: print "reading attribute '%s'" %key
1022           
1023        attributes[key] = array(fid.variables[key])
1024   
1025   
1026    try:
1027        geo_reference = Geo_reference(NetCDFObject=fid)
1028    except AttributeError, e:
1029        geo_reference = None
1030   
1031    fid.close()
1032   
1033    return pointlist, attributes, geo_reference
1034
1035
1036def _read_csv_file(file_name, verbose=False):
1037    """Read .csv file
1038   
1039    Return a dic of array of points, and dic of array of attribute
1040    eg
1041    dic['points'] = [[1.0,2.0],[3.0,5.0]]
1042    dic['attributelist']['elevation'] = [[7.0,5.0]
1043    """
1044   
1045    file_pointer = open(file_name)
1046    header, file_pointer = _read_csv_file_header(file_pointer)
1047    try:
1048        pointlist, att_dict, geo_ref, file_pointer = \
1049                   _read_csv_file_blocking( \
1050                file_pointer,
1051                header,
1052                max_read_lines=1e30) #If the file is bigger that this, block..
1053    except ANUGAError:
1054        file_pointer.close()
1055        raise 
1056    file_pointer.close()
1057    return pointlist, att_dict, geo_ref   
1058
1059CSV_DELIMITER = ','
1060def _read_csv_file_header(file_pointer,
1061                          delimiter=CSV_DELIMITER,
1062                          verbose=False):
1063
1064    """Read the header of a .csv file
1065    Return a list of the header names
1066    """
1067    line = file_pointer.readline()
1068    header = clean_line(line, delimiter)
1069    return header, file_pointer
1070
1071def _read_csv_file_blocking(file_pointer, header,
1072                            delimiter=CSV_DELIMITER,
1073                            max_read_lines=MAX_READ_LINES,
1074                            verbose=False):
1075   
1076
1077    """
1078    Read the body of a .csv file.
1079    header: The list header of the csv file, with the x and y labels.
1080    """
1081    points = []
1082    pointattributes = []
1083    att_dict = {}
1084
1085    #This is to remove the x and y headers.
1086    header = header[:]
1087    try:
1088        x_header = header.pop(0)
1089        y_header = header.pop(0)
1090    except IndexError:
1091        # if there are not two columns this will occur.
1092        # eg if it is a space seperated file
1093        raise SyntaxError
1094   
1095    read_lines = 0
1096    while read_lines<max_read_lines:
1097        line = file_pointer.readline()
1098        #print "line",line
1099        numbers = clean_line(line,delimiter)
1100        if len(numbers) <= 1:
1101            break
1102        if line[0] == '#':
1103            continue
1104        read_lines += 1
1105        try:
1106            x = float(numbers[0])
1107            y = float(numbers[1])
1108            points.append([x,y])
1109            numbers.pop(0)
1110            numbers.pop(0)
1111            if len(header) != len(numbers):
1112                file_pointer.close() 
1113                msg = "File load error.  There might be a problem with the file header"
1114                raise SyntaxError, msg
1115            for i,num in enumerate(numbers):
1116                num.strip()
1117                if num != '\n' and num != '':
1118                    #attributes.append(float(num))
1119                    att_dict.setdefault(header[i],[]).append(float(num))
1120        #except IOError:           
1121        except ValueError:
1122            raise SyntaxError
1123    if points == []:
1124        raise StopIteration
1125       
1126   
1127    pointlist = array(points).astype(Float)
1128    for key in att_dict.keys():
1129        att_dict[key] = array(att_dict[key]).astype(Float)
1130
1131    #Do stuff here so the info is in lat's and longs
1132    geo_ref = None
1133    x_header = lower(x_header[:3])
1134    y_header = lower(y_header[:3])
1135    if (x_header == 'lon' or  x_header == 'lat') and \
1136       (y_header == 'lon' or  y_header == 'lat'):
1137        if x_header == 'lon':
1138            longitudes = ravel(pointlist[:,0:1])
1139            latitudes = ravel(pointlist[:,1:])
1140        else:
1141            latitudes = ravel(pointlist[:,0:1])
1142            longitudes = ravel(pointlist[:,1:])
1143       
1144        pointlist, geo_ref = _set_using_lat_long(latitudes,
1145                                                 longitudes,
1146                                                 geo_reference=None,
1147                                                 data_points=None,
1148                                                 points_are_lats_longs=False)
1149    return pointlist, att_dict, geo_ref, file_pointer
1150
1151def _read_pts_file_header(fid, verbose=False):
1152
1153    """Read the geo_reference and number_of_points from a .pts file
1154    """
1155   
1156    keys = fid.variables.keys()
1157    try:
1158        keys.remove('points')
1159    except IOError, e:       
1160        fid.close()   
1161        msg = 'Expected keyword "points" but could not find it'
1162        raise IOError, msg
1163    if verbose: print 'Got %d variables: %s' %(len(keys), keys)
1164   
1165    try:
1166        geo_reference = Geo_reference(NetCDFObject=fid)
1167    except AttributeError, e:
1168        geo_reference = None
1169
1170    return geo_reference, keys, fid.dimensions['number_of_points']
1171
1172def _read_pts_file_blocking(fid, start_row, fin_row, keys):
1173                            #verbose=False):
1174   
1175
1176    """
1177    Read the body of a .csv file.
1178    header: The list header of the csv file, with the x and y labels.
1179    """
1180   
1181    pointlist = array(fid.variables['points'][start_row:fin_row])
1182   
1183    attributes = {}
1184    for key in keys:
1185        attributes[key] = array(fid.variables[key][start_row:fin_row])
1186
1187    return pointlist, attributes
1188   
1189   
1190def _read_xya_file(fd, delimiter):
1191    points = []
1192    pointattributes = []
1193    title = fd.readline()
1194    att_names = clean_line(title,delimiter)
1195    att_dict = {}
1196    line = fd.readline()
1197    numbers = clean_line(line,delimiter)
1198   
1199    while len(numbers) > 1 and line[0] <> '#':
1200        if numbers != []:
1201            try:
1202                x = float(numbers[0])
1203                y = float(numbers[1])
1204                points.append([x,y])
1205                numbers.pop(0)
1206                numbers.pop(0)
1207                if len(att_names) != len(numbers):
1208                    fd.close()
1209                    # It might not be a problem with the title
1210                    #raise TitleAmountError
1211                    raise IOError
1212                for i,num in enumerate(numbers):
1213                    num.strip()
1214                    if num != '\n' and num != '':
1215                        #attributes.append(float(num))
1216                        att_dict.setdefault(att_names[i],[]).append(float(num))
1217            except ValueError:
1218                raise SyntaxError
1219        line = fd.readline()
1220        numbers = clean_line(line,delimiter)
1221   
1222    if line == '':
1223        geo_reference = None
1224    else:
1225        geo_reference = Geo_reference(ASCIIFile=fd,read_title=line)
1226       
1227   
1228    pointlist = array(points).astype(Float)
1229    for key in att_dict.keys():
1230        att_dict[key] = array(att_dict[key]).astype(Float)
1231   
1232    return pointlist, att_dict, geo_reference
1233
1234def _write_pts_file(file_name,
1235                    write_data_points,
1236                    write_attributes=None, 
1237                    write_geo_reference=None):
1238    """
1239    Write .pts NetCDF file   
1240
1241    NOTE: Below might not be valid ask Duncan : NB 5/2006
1242   
1243    WARNING: This function mangles the point_atts data structure
1244    #F??ME: (DSG)This format has issues.
1245    # There can't be an attribute called points
1246    # consider format change
1247    # method changed by NB not sure if above statement is correct
1248   
1249    should create new test for this
1250    legal_keys = ['pointlist', 'attributelist', 'geo_reference']
1251    for key in point_atts.keys():
1252        msg = 'Key %s is illegal. Valid keys are %s' %(key, legal_keys)
1253        assert key in legal_keys, msg
1254    """   
1255    from Scientific.IO.NetCDF import NetCDFFile
1256    # NetCDF file definition
1257    outfile = NetCDFFile(file_name, 'w')
1258   
1259    #Create new file
1260    outfile.institution = 'Geoscience Australia'
1261    outfile.description = 'NetCDF format for compact and portable storage ' +\
1262                          'of spatial point data'
1263   
1264    # dimension definitions
1265    shape = write_data_points.shape[0]
1266    outfile.createDimension('number_of_points', shape) 
1267    outfile.createDimension('number_of_dimensions', 2) #This is 2d data
1268   
1269    # variable definition
1270    outfile.createVariable('points', Float, ('number_of_points',
1271                                             'number_of_dimensions'))
1272
1273    #create variables 
1274    outfile.variables['points'][:] = write_data_points #.astype(Float32)
1275
1276    if write_attributes is not None:
1277        for key in write_attributes.keys():
1278            outfile.createVariable(key, Float, ('number_of_points',))
1279            outfile.variables[key][:] = write_attributes[key] #.astype(Float32)
1280       
1281    if write_geo_reference is not None:
1282        write_NetCDF_georeference(write_geo_reference, outfile)
1283       
1284    outfile.close() 
1285 
1286
1287
1288def _write_xya_file(file_name,
1289                    write_data_points,
1290                    write_attributes=None, 
1291                    write_geo_reference=None, 
1292                    delimiter=','):
1293    """
1294    export a file, file_name, with the xya format
1295   
1296    """
1297    points = write_data_points
1298    pointattributes = write_attributes
1299   
1300    fd = open(file_name,'w')
1301    titlelist = ""
1302    if pointattributes is not None:   
1303        for title in pointattributes.keys():
1304            titlelist = titlelist + title + delimiter
1305        titlelist = titlelist[0:-len(delimiter)] # remove the last delimiter
1306    fd.write(titlelist+"\n")
1307   
1308    #<vertex #> <x> <y> [attributes]
1309    for i, vert in enumerate( points):
1310
1311
1312        if pointattributes is not None:           
1313            attlist = ","
1314            for att in pointattributes.keys():
1315                attlist = attlist + str(pointattributes[att][i])+ delimiter
1316            attlist = attlist[0:-len(delimiter)] # remove the last delimiter
1317            attlist.strip()
1318        else:
1319            attlist = ''
1320
1321        fd.write(str(vert[0]) + delimiter +
1322                 str(vert[1]) + attlist + "\n")
1323
1324    if  write_geo_reference is not None:
1325        write_geo_reference = ensure_geo_reference(write_geo_reference)
1326        write_geo_reference.write_ASCII(fd)
1327    fd.close()
1328
1329
1330def _write_csv_file(file_name,
1331                    write_data_points,
1332                    write_attributes=None,
1333                    as_lat_long=False,
1334                    delimiter=','):
1335    """
1336    export a file, file_name, with the csv format
1337   
1338    """
1339    points = write_data_points
1340    pointattributes = write_attributes
1341   
1342    fd = open(file_name,'w')
1343    if as_lat_long:
1344        titlelist = "latitude" + delimiter + "longitude"  + delimiter
1345    else:
1346        titlelist = "x" + delimiter + "y"  + delimiter
1347    if pointattributes is not None:   
1348        for title in pointattributes.keys():
1349            titlelist = titlelist + title + delimiter
1350        titlelist = titlelist[0:-len(delimiter)] # remove the last delimiter
1351    fd.write(titlelist+"\n")
1352   
1353    # <x/lat> <y/long> [attributes]
1354    for i, vert in enumerate( points):
1355
1356
1357        if pointattributes is not None:           
1358            attlist = ","
1359            for att in pointattributes.keys():
1360                attlist = attlist + str(pointattributes[att][i])+ delimiter
1361            attlist = attlist[0:-len(delimiter)] # remove the last delimiter
1362            attlist.strip()
1363        else:
1364            attlist = ''
1365
1366        fd.write(str(vert[0]) + delimiter +
1367                 str(vert[1]) + attlist + "\n")
1368
1369    fd.close()
1370 
1371def _write_urs_file(file_name,
1372                    points,
1373                    delimiter=' '):
1374    """
1375    export a file, file_name, with the urs format
1376    the data points are in lats and longs
1377   
1378    """   
1379    fd = open(file_name,'w')   
1380    fd.write(str(len(points))+"\n")   
1381    # <lat> <long> <id#>
1382    for i, vert in enumerate( points):
1383        fd.write(str(round(vert[0],7)) + delimiter + \
1384                 str(round(vert[1],7)) + delimiter +str(i)+ "\n")
1385    fd.close()
1386   
1387def _point_atts2array(point_atts):
1388    point_atts['pointlist'] = array(point_atts['pointlist']).astype(Float)
1389   
1390    for key in point_atts['attributelist'].keys():
1391        point_atts['attributelist'][key]=\
1392                array(point_atts['attributelist'][key]).astype(Float)
1393    return point_atts
1394
1395   
1396
1397
1398def geospatial_data2points_dictionary(geospatial_data):
1399    """Convert geospatial data to points_dictionary
1400    """
1401
1402    points_dictionary = {}
1403    points_dictionary['pointlist'] = geospatial_data.data_points
1404
1405    points_dictionary['attributelist'] = {}
1406
1407    for attribute_name in geospatial_data.attributes.keys():
1408        val = geospatial_data.attributes[attribute_name]
1409        points_dictionary['attributelist'][attribute_name] = val
1410
1411    points_dictionary['geo_reference'] = geospatial_data.geo_reference
1412
1413    return points_dictionary
1414
1415   
1416def points_dictionary2geospatial_data(points_dictionary):
1417    """Convert points_dictionary to geospatial data object
1418    """
1419
1420    msg = 'Points dictionary must have key pointlist' 
1421    assert points_dictionary.has_key('pointlist'), msg
1422
1423    msg = 'Points dictionary must have key attributelist'     
1424    assert points_dictionary.has_key('attributelist'), msg       
1425
1426    if points_dictionary.has_key('geo_reference'):
1427        geo = points_dictionary['geo_reference']
1428    else:
1429        geo = None
1430   
1431    return Geospatial_data(points_dictionary['pointlist'],
1432                           points_dictionary['attributelist'],
1433                           geo_reference = geo)
1434
1435def clean_line(line,delimiter):     
1436    """Remove whitespace
1437    """
1438    #print ">%s" %line
1439    line = line.strip()
1440    #print "stripped>%s" %line
1441    numbers = line.split(delimiter)
1442    i = len(numbers) - 1
1443    while i >= 0:
1444        if numbers[i] == '':
1445            numbers.pop(i)
1446        else:
1447            numbers[i] = numbers[i].strip()
1448       
1449        i += -1
1450    #for num in numbers:
1451    #    print "num>%s<" %num
1452    return numbers
1453           
1454def ensure_absolute(points, geo_reference=None):
1455    """
1456    This function inputs several formats and
1457    outputs one format. - a numeric array of absolute points.
1458
1459    Inputed formats are;
1460    points: List or numeric array of coordinate pairs [xi, eta] of
1461              points or geospatial object or points file name
1462
1463    mesh_origin: A geo_reference object or 3-tuples consisting of
1464                 UTM zone, easting and northing.
1465                 If specified vertex coordinates are assumed to be
1466                 relative to their respective origins.
1467    """
1468    if isinstance(points,type('')):
1469        #It's a string
1470        #assume it is a point file
1471        points = Geospatial_data(file_name = points)
1472       
1473    if isinstance(points,Geospatial_data):
1474        points = points.get_data_points( \
1475                absolute = True)
1476        msg = "Use a Geospatial_data object or a mesh origin. Not both."
1477        assert geo_reference == None, msg
1478           
1479    else:
1480        points = ensure_numeric(points, Float)
1481    if geo_reference is None:
1482        geo = None #Geo_reference()
1483    else:
1484        if isinstance(geo_reference, Geo_reference):
1485            geo = geo_reference
1486        else:
1487            geo = Geo_reference(geo_reference[0],
1488                                geo_reference[1],
1489                                geo_reference[2])
1490        points = geo.get_absolute(points)
1491    return points
1492     
1493
1494def ensure_geospatial(points, geo_reference=None):
1495    """
1496    This function inputs several formats and
1497    outputs one format. - a geospatial_data instance.
1498
1499    Inputed formats are;
1500    points: List or numeric array of coordinate pairs [xi, eta] of
1501              points or geospatial object
1502
1503    mesh_origin: A geo_reference object or 3-tuples consisting of
1504                 UTM zone, easting and northing.
1505                 If specified vertex coordinates are assumed to be
1506                 relative to their respective origins.
1507    """
1508    if isinstance(points,Geospatial_data):
1509        msg = "Use a Geospatial_data object or a mesh origin. Not both."
1510        assert geo_reference == None, msg
1511        return points   
1512    else:
1513        points = ensure_numeric(points, Float)
1514    if geo_reference is None:
1515        geo = None
1516    else:
1517        if isinstance(geo_reference, Geo_reference):
1518            geo = geo_reference
1519        else:
1520            geo = Geo_reference(geo_reference[0],
1521                                geo_reference[1],
1522                                geo_reference[2])
1523    points = Geospatial_data(data_points=points, geo_reference=geo)       
1524    return points
1525
1526#def file2xya(filename):
1527   
1528#    G = Geospatial_data(filename)
1529#    G.export_points_file()
1530
1531             
1532
1533         
1534if __name__ == "__main__":
1535    pass
1536   
Note: See TracBrowser for help on using the repository browser.