source: branches/numpy/anuga/coordinate_transforms/geo_reference.py @ 6883

Last change on this file since 6883 was 6689, checked in by rwilson, 16 years ago

Back-merge from Numeric trunk to numpy branch.

File size: 17.4 KB
RevLine 
[2488]1"""
[4451]2
[2488]3"""
4
5
6#FIXME: Ensure that all attributes of a georef are treated everywhere
7#and unit test
8
[2593]9import types, sys
[6441]10import copy
11
[3514]12from anuga.utilities.numerical_tools import ensure_numeric
[6360]13from anuga.utilities.anuga_exceptions import ANUGAError, TitleError, \
14                                             ParsingError, ShapeError
[6304]15from anuga.config import netcdf_float, netcdf_int, netcdf_float32
[2488]16
[6304]17import numpy as num
[2941]18
[6149]19
[2602]20DEFAULT_ZONE = -1
[6360]21TITLE = '#geo reference' + "\n" # this title is referred to in the test format
[2488]22
[4387]23DEFAULT_PROJECTION = 'UTM'
24DEFAULT_DATUM = 'wgs84'
25DEFAULT_UNITS = 'm'
26DEFAULT_FALSE_EASTING = 500000
[6360]27DEFAULT_FALSE_NORTHING = 10000000    # Default for southern hemisphere
[4387]28
[6360]29
30##
31# @brief A class for ...
[2488]32class Geo_reference:
33    """
[6360]34    Attributes of the Geo_reference class:
35        .zone           The UTM zone (default is -1)
36        .false_easting  ??
37        .false_northing ??
38        .datum          The Datum used (default is wgs84)
39        .projection     The projection used (default is 'UTM')
40        .units          The units of measure used (default metres)
41        .xllcorner      The X coord of origin (default is 0.0 wrt UTM grid)
42        .yllcorner      The y coord of origin (default is 0.0 wrt UTM grid)
43        .is_absolute    ??
44
[2488]45    """
46
[6360]47    ##
48    # @brief Instantiate an instance of class Geo_reference.
49    # @param zone The UTM zone.
50    # @param xllcorner X coord of origin of georef.
51    # @param yllcorner Y coord of origin of georef.
52    # @param datum ??
53    # @param projection The projection used (default UTM).
54    # @param units Units used in measuring distance (default m).
55    # @param false_easting ??
56    # @param false_northing ??
57    # @param NetCDFObject NetCDF file *handle* to write to.
58    # @param ASCIIFile ASCII text file *handle* to write to.
59    # @param read_title Title of the georeference text.
[2488]60    def __init__(self,
[6360]61                 zone=DEFAULT_ZONE,
62                 xllcorner=0.0,
63                 yllcorner=0.0,
64                 datum=DEFAULT_DATUM,
65                 projection=DEFAULT_PROJECTION,
66                 units=DEFAULT_UNITS,
67                 false_easting=DEFAULT_FALSE_EASTING,
68                 false_northing=DEFAULT_FALSE_NORTHING,
[2488]69                 NetCDFObject=None,
70                 ASCIIFile=None,
71                 read_title=None):
72        """
73        input:
[6360]74        NetCDFObject - a handle to the netCDF file to be written to
[2488]75        ASCIIFile - a handle to the text file
76        read_title - the title of the georeference text, if it was read in.
[2941]77         If the function that calls this has already read the title line,
78         it can't unread it, so this info has to be passed.
79         If you know of a way to unread this info, then tell us.
[2488]80
[4387]81         Note, the text file only saves a sub set of the info the
82         points file does.  Currently the info not written in text
83         must be the default info, since ANUGA assumes it isn't
84         changing.
[6360]85        """
86
[4451]87        if zone is None:
88            zone = DEFAULT_ZONE
[2488]89        self.false_easting = false_easting
[6360]90        self.false_northing = false_northing
[2488]91        self.datum = datum
92        self.projection = projection
[6360]93        self.zone = zone
[2488]94        self.units = units
95        self.xllcorner = xllcorner
[6553]96        self.yllcorner = yllcorner       
97           
[2488]98        if NetCDFObject is not None:
99            self.read_NetCDF(NetCDFObject)
[6360]100
[2488]101        if ASCIIFile is not None:
102            self.read_ASCII(ASCIIFile, read_title=read_title)
103
[6553]104           
105        # Set flag for absolute points (used by get_absolute)   
106        self.absolute = num.allclose([self.xllcorner, self.yllcorner], 0)
107           
108
[2488]109    def get_xllcorner(self):
110        return self.xllcorner
[6360]111
112    ##
113    # @brief Get the Y coordinate of the origin of this georef.
[2488]114    def get_yllcorner(self):
115        return self.yllcorner
[6360]116
117    ##
118    # @brief Get the zone of this georef.
[2488]119    def get_zone(self):
120        return self.zone
[6360]121
122    ##
123    # @brief Write <something> to an open NetCDF file.
124    # @param outfile Handle to open NetCDF file.
[2488]125    def write_NetCDF(self, outfile):
[6360]126        outfile.xllcorner = self.xllcorner
[2488]127        outfile.yllcorner = self.yllcorner
128        outfile.zone = self.zone
129
130        outfile.false_easting = self.false_easting
131        outfile.false_northing = self.false_northing
132
[6360]133        outfile.datum = self.datum
[2488]134        outfile.projection = self.projection
135        outfile.units = self.units
136
[6360]137    ##
138    # @brief Read data from an open NetCDF file.
139    # @param infile Handle to open NetCDF file.
[2488]140    def read_NetCDF(self, infile):
141        self.xllcorner = infile.xllcorner[0]
[6360]142        self.yllcorner = infile.yllcorner[0]
[2488]143        self.zone = infile.zone[0]
144
145        # Fix some assertion failures
[6304]146        if isinstance(self.zone, num.ndarray) and self.zone.shape == ():
[2488]147            self.zone = self.zone[0]
[6360]148        if (isinstance(self.xllcorner, num.ndarray) and
149                self.xllcorner.shape == ()):
[2488]150            self.xllcorner = self.xllcorner[0]
[6360]151        if (isinstance(self.yllcorner, num.ndarray) and
152                self.yllcorner.shape == ()):
[2488]153            self.yllcorner = self.yllcorner[0]
154
[6304]155        assert (self.xllcorner.dtype.kind in num.typecodes['Float'] or
156                self.xllcorner.dtype.kind in num.typecodes['Integer'])
157        assert (self.yllcorner.dtype.kind in num.typecodes['Float'] or
158                self.yllcorner.dtype.kind in num.typecodes['Integer'])
159        assert (self.zone.dtype.kind in num.typecodes['Integer'])
[6360]160
[4389]161        try:
162            self.false_easting = infile.false_easting[0]
163            self.false_northing = infile.false_northing[0]
[6360]164
165            self.datum = infile.datum
[4389]166            self.projection = infile.projection
167            self.units = infile.units
168        except:
169            pass
[6360]170
171        if self.false_easting != DEFAULT_FALSE_EASTING:
172            print "WARNING: False easting of %f specified." % self.false_easting
173            print "Default false easting is %f." % DEFAULT_FALSE_EASTING
[4451]174            print "ANUGA does not correct for differences in False Eastings."
[6360]175
176        if self.false_northing != DEFAULT_FALSE_NORTHING:
177            print ("WARNING: False northing of %f specified."
178                   % self.false_northing)
179            print "Default false northing is %f." % DEFAULT_FALSE_NORTHING
[4451]180            print "ANUGA does not correct for differences in False Northings."
[6360]181
182        if self.datum.upper() != DEFAULT_DATUM.upper():
183            print "WARNING: Datum of %s specified." % self.datum
184            print "Default Datum is %s." % DEFAULT_DATUM
[4451]185            print "ANUGA does not correct for differences in datums."
[6360]186
187        if self.projection.upper() != DEFAULT_PROJECTION.upper():
188            print "WARNING: Projection of %s specified." % self.projection
189            print "Default Projection is %s." % DEFAULT_PROJECTION
[4451]190            print "ANUGA does not correct for differences in Projection."
[6360]191
192        if self.units.upper() != DEFAULT_UNITS.upper():
193            print "WARNING: Units of %s specified." % self.units
194            print "Default units is %s." % DEFAULT_UNITS
[4451]195            print "ANUGA does not correct for differences in units."
[6360]196
197################################################################################
198# ASCII files with geo-refs are currently not used
199################################################################################
200
201    ##
202    # @brief Write georef data to an open text file.
203    # @param fd Handle to open text file.
[2488]204    def write_ASCII(self, fd):
205        fd.write(TITLE)
[6360]206        fd.write(str(self.zone) + "\n")
207        fd.write(str(self.xllcorner) + "\n")
[2488]208        fd.write(str(self.yllcorner) + "\n")
209
[6360]210    ##
211    # @brief Read georef data from an open text file.
212    # @param fd Handle to open text file.
213    def read_ASCII(self, fd, read_title=None):
[2942]214        try:
215            if read_title == None:
[6360]216                read_title = fd.readline()     # remove the title line
[2942]217            if read_title[0:2].upper() != TITLE[0:2].upper():
[6360]218                msg = ('File error.  Expecting line: %s.  Got this line: %s'
219                       % (TITLE, read_title))
220                raise TitleError, msg
[2942]221            self.zone = int(fd.readline())
222            self.xllcorner = float(fd.readline())
223            self.yllcorner = float(fd.readline())
224        except SyntaxError:
[6360]225            msg = 'File error.  Got syntax error while parsing geo reference'
226            raise ParsingError, msg
227
[2488]228        # Fix some assertion failures
[6304]229        if isinstance(self.zone, num.ndarray) and self.zone.shape == ():
[2488]230            self.zone = self.zone[0]
[6360]231        if (isinstance(self.xllcorner, num.ndarray) and
232                self.xllcorner.shape == ()):
[2488]233            self.xllcorner = self.xllcorner[0]
[6360]234        if (isinstance(self.yllcorner, num.ndarray) and
235                self.yllcorner.shape == ()):
[2488]236            self.yllcorner = self.yllcorner[0]
[6304]237
[6553]238        assert (type(self.xllcorner) == types.FloatType)
239        assert (type(self.yllcorner) == types.FloatType)
240        assert (type(self.zone) == types.IntType)
241
[6360]242################################################################################
[6304]243
[6428]244    ##
245    # @brief Change points to be absolute wrt new georef 'points_geo_ref'.
246    # @param points The points to change.
247    # @param points_geo_ref The new georef to make points absolute wrt.
248    # @return The changed points.
249    # @note If 'points' is a list then a changed list is returned.
[6360]250    def change_points_geo_ref(self, points, points_geo_ref=None):
[6553]251        """Change the geo reference of a list or Numeric array of points to
[2488]252        be this reference.(The reference used for this object)
[6553]253        If the points do not have a geo ref, assume 'absolute' values
[2488]254        """
[6553]255        import copy
256       
[6428]257        # remember if we got a list
258        is_list = isinstance(points, list)
[2594]259
[6553]260        points = ensure_numeric(points, num.float)
[6360]261
[6553]262        # sanity checks
[2594]263        if len(points.shape) == 1:
[6553]264            #One point has been passed
[2594]265            msg = 'Single point must have two elements'
266            assert len(points) == 2, msg
[6149]267            points = num.reshape(points, (1,2))
[2594]268
[6553]269        msg = 'Points array must be two dimensional.\n'
270        msg += 'I got %d dimensions' %len(points.shape)
[5665]271        assert len(points.shape) == 2, msg
272
[6553]273        msg = 'Input must be an N x 2 array or list of (x,y) values. '
274        msg += 'I got an %d x %d array' %points.shape   
275        assert points.shape[1] == 2, msg               
[3071]276
[6553]277        # FIXME (Ole): Could also check if zone, xllcorner, yllcorner
278        # are identical in the two geo refs.   
[2594]279        if points_geo_ref is not self:
[5736]280            # If georeferences are different
[6553]281            points = copy.copy(points) # Don't destroy input                   
[2594]282            if not points_geo_ref is None:
[5736]283                # Convert points to absolute coordinates
[6553]284                points[:,0] += points_geo_ref.xllcorner
285                points[:,1] += points_geo_ref.yllcorner
286       
[5736]287            # Make points relative to primary geo reference
[6553]288            points[:,0] -= self.xllcorner
[2594]289            points[:,1] -= self.yllcorner
290
[2488]291        if is_list:
292            points = points.tolist()
[6553]293           
[2488]294        return points
[3129]295
296    def is_absolute(self):
[6553]297        """Return True if xllcorner==yllcorner==0 indicating that points
298        in question are absolute.
299        """
300       
301        # FIXME(Ole): It is unfortunate that decision about whether points
302        # are absolute or not lies with the georeference object. Ross pointed this out.
303        # Moreover, this little function is responsible for a large fraction of the time
304        # using in data fitting (something in like 40 - 50%.
305        # This was due to the repeated calls to allclose.
306        # With the flag method fitting is much faster (18 Mar 2009).
[3129]307
[6553]308        # FIXME(Ole): HACK to be able to reuse data already cached (18 Mar 2009).
309        # Remove at some point
310        if not hasattr(self, 'absolute'):
311            self.absolute = num.allclose([self.xllcorner, self.yllcorner], 0)
312           
313        # Return absolute flag   
314        return self.absolute
[3129]315
[2488]316    def get_absolute(self, points):
[6553]317        """Given a set of points geo referenced to this instance,
[2488]318        return the points as absolute values.
319        """
[3071]320
[6553]321        is_list = False
322        if type(points) == types.ListType:
323            is_list = True
[3071]324
[6553]325        points = ensure_numeric(points, num.float)
[3071]326        if len(points.shape) == 1:
[6553]327            # One point has been passed
328            msg = 'Single point must have two elements'
[3134]329            if not len(points) == 2:
[6553]330                raise ShapeError, msg   
[3071]331
[6553]332
333        msg = 'Input must be an N x 2 array or list of (x,y) values. '
334        msg += 'I got an %d x %d array' %points.shape   
[3134]335        if not points.shape[1] == 2:
[6553]336            raise ShapeError, msg   
337           
338       
[3134]339        # Add geo ref to points
340        if not self.is_absolute():
[6553]341            points = copy.copy(points) # Don't destroy input                   
342            points[:,0] += self.xllcorner
[3134]343            points[:,1] += self.yllcorner
[6360]344
[6553]345       
[2488]346        if is_list:
347            points = points.tolist()
[6553]348             
[2488]349        return points
350
[6360]351    ##
352    # @brief Convert points to relative measurement.
353    # @param points Points to convert to relative measurements.
354    # @return A set of points relative to the geo_reference instance.
[5288]355    def get_relative(self, points):
[6360]356        """Given a set of points in absolute UTM coordinates,
[5288]357        make them relative to this geo_reference instance,
358        return the points as relative values.
359
360        This is the inverse of get_absolute.
361        """
362
[6442]363        # remember if we got a list
364        is_list = isinstance(points, list)
[5288]365
[6553]366        points = ensure_numeric(points, num.float)
[5288]367        if len(points.shape) == 1:
368            #One point has been passed
[6553]369            msg = 'Single point must have two elements'
[5288]370            if not len(points) == 2:
[6553]371                raise ShapeError, msg   
[5288]372
[6360]373        if not points.shape[1] == 2:
374            msg = ('Input must be an N x 2 array or list of (x,y) values. '
375                   'I got an %d x %d array' % points.shape)
[6553]376            raise ShapeError, msg   
[5288]377
378        # Subtract geo ref from points
379        if not self.is_absolute():
[6553]380            points = copy.copy(points) # Don't destroy input                           
381            points[:,0] -= self.xllcorner
[5288]382            points[:,1] -= self.yllcorner
383
384        if is_list:
385            points = points.tolist()
[6553]386             
[5288]387        return points
388
[6360]389    ##
390    # @brief ??
391    # @param other ??
[2594]392    def reconcile_zones(self, other):
[4180]393        if other is None:
394            other = Geo_reference()
[6360]395        if (self.zone == other.zone or
396            self.zone == DEFAULT_ZONE and
397            other.zone == DEFAULT_ZONE):
[2594]398            pass
399        elif self.zone == DEFAULT_ZONE:
400            self.zone = other.zone
401        elif other.zone == DEFAULT_ZONE:
[6360]402            other.zone = self.zone
403        else:
404            msg = ('Geospatial data must be in the same '
405                   'ZONE to allow reconciliation. I got zone %d and %d'
406                   % (self.zone, other.zone))
[2608]407            raise ANUGAError, msg
[6360]408
[2488]409    #def easting_northing2geo_reffed_point(self, x, y):
410    #    return [x-self.xllcorner, y - self.xllcorner]
411
412    #def easting_northing2geo_reffed_points(self, x, y):
413    #    return [x-self.xllcorner, y - self.xllcorner]
414
[6360]415    ##
416    # @brief Get origin of this geo_reference.
417    # @return (zone, xllcorner, yllcorner).
[2488]418    def get_origin(self):
419        return (self.zone, self.xllcorner, self.yllcorner)
[6360]420
421    ##
422    # @brief Get a string representation of this geo_reference instance.
[2488]423    def __repr__(self):
[6360]424        return ('(zone=%i easting=%f, northing=%f)'
425                % (self.zone, self.xllcorner, self.yllcorner))
426
427    ##
428    # @brief Compare two geo_reference instances.
429    # @param self This geo_reference instance.
430    # @param other Another geo_reference instance to compare against.
431    # @return 0 if instances have the same attributes, else 1.
432    # @note Attributes are: zone, xllcorner, yllcorner.
433    def __cmp__(self, other):
[2488]434        # FIXME (DSG) add a tolerence
435        if other is None:
436            return 1
437        cmp = 0
438        if not (self.xllcorner == self.xllcorner):
439            cmp = 1
440        if not (self.yllcorner == self.yllcorner):
441            cmp = 1
442        if not (self.zone == self.zone):
443            cmp = 1
444        return cmp
[4387]445
[6360]446
447##
448# @brief Write a geo_reference to a NetCDF file (usually SWW).
449# @param origin A georef instance or parameters to create a georef instance.
450# @param outfile Path to file to write.
451# @return A normalized geo_reference.
[4387]452def write_NetCDF_georeference(origin, outfile):
[6360]453    """Write georeference info to a netcdf file, usually sww.
[4387]454
[6360]455    The origin can be a georef instance or parameters for a geo_ref instance
[4387]456
457    outfile is the name of the file to be written to.
458    """
[6360]459
[4451]460    geo_ref = ensure_geo_reference(origin)
461    geo_ref.write_NetCDF(outfile)
462    return geo_ref
463
[6360]464
465##
466# @brief Convert an object to a georeference instance.
467# @param origin A georef instance or (zone, xllcorner, yllcorner)
468# @return A georef object, or None if 'origin' was None.
[4451]469def ensure_geo_reference(origin):
470    """
471    Given a list/tuple of zone, xllcorner and yllcorner of a geo-ref object,
472    return a geo ref object.
473
474    If the origin is None, return None, so calling this function doesn't
475    effect code logic
476    """
[6360]477
478    if isinstance(origin, Geo_reference):
[4387]479        geo_ref = origin
[4451]480    elif origin is None:
481        geo_ref = None
[4387]482    else:
[6360]483        geo_ref = apply(Geo_reference, origin)
484
[4387]485    return geo_ref
[4451]486
[6360]487
[2488]488#-----------------------------------------------------------------------
489
490if __name__ == "__main__":
491    pass
Note: See TracBrowser for help on using the repository browser.