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

"""

"""


#FIXME: Ensure that all attributes of a georef are treated everywhere
#and unit test

import types, sys
from anuga.utilities.numerical_tools import ensure_numeric
from anuga.utilities.anuga_exceptions import ANUGAError, TitleError, \
                                             ParsingError, ShapeError
from anuga.config import netcdf_float, netcdf_int, netcdf_float32

import numpy as num


DEFAULT_ZONE = -1
TITLE = '#geo reference' + "\n" # this title is referred to in the test format

DEFAULT_PROJECTION = 'UTM'
DEFAULT_DATUM = 'wgs84'
DEFAULT_UNITS = 'm'
DEFAULT_FALSE_EASTING = 500000
DEFAULT_FALSE_NORTHING = 10000000    # Default for southern hemisphere


##
# @brief A class for ...
class Geo_reference:
    """
    Attributes of the Geo_reference class:
        .zone           The UTM zone (default is -1)
        .false_easting  ??
        .false_northing ??
        .datum          The Datum used (default is wgs84)
        .projection     The projection used (default is 'UTM')
        .units          The units of measure used (default metres)
        .xllcorner      The X coord of origin (default is 0.0 wrt UTM grid)
        .yllcorner      The y coord of origin (default is 0.0 wrt UTM grid)
        .is_absolute    ??

    """

    ##
    # @brief Instantiate an instance of class Geo_reference.
    # @param zone The UTM zone.
    # @param xllcorner X coord of origin of georef.
    # @param yllcorner Y coord of origin of georef.
    # @param datum ??
    # @param projection The projection used (default UTM).
    # @param units Units used in measuring distance (default m).
    # @param false_easting ??
    # @param false_northing ??
    # @param NetCDFObject NetCDF file *handle* to write to.
    # @param ASCIIFile ASCII text file *handle* to write to.
    # @param read_title Title of the georeference text.
    def __init__(self,
                 zone=DEFAULT_ZONE,
                 xllcorner=0.0,
                 yllcorner=0.0,
                 datum=DEFAULT_DATUM,
                 projection=DEFAULT_PROJECTION,
                 units=DEFAULT_UNITS,
                 false_easting=DEFAULT_FALSE_EASTING,
                 false_northing=DEFAULT_FALSE_NORTHING,
                 NetCDFObject=None,
                 ASCIIFile=None,
                 read_title=None):
        """
        input:
        NetCDFObject - a handle to the netCDF file to be written to
        ASCIIFile - a handle to the text file
        read_title - the title of the georeference text, if it was read in.
         If the function that calls this has already read the title line,
         it can't unread it, so this info has to be passed.
         If you know of a way to unread this info, then tell us.

         Note, the text file only saves a sub set of the info the
         points file does.  Currently the info not written in text
         must be the default info, since ANUGA assumes it isn't
         changing.
        """

        if zone is None:
            zone = DEFAULT_ZONE
        self.false_easting = false_easting
        self.false_northing = false_northing
        self.datum = datum
        self.projection = projection
        self.zone = zone
        self.units = units
        self.xllcorner = xllcorner
        self.yllcorner = yllcorner
        #self.is_absolute = num.allclose([self.xllcorner, self.yllcorner], 0)

        if NetCDFObject is not None:
            self.read_NetCDF(NetCDFObject)

        if ASCIIFile is not None:
            self.read_ASCII(ASCIIFile, read_title=read_title)

#    # Might be better to have this method instead of the following 3.
#    def get_origin(self):
#        return (self.zone, self.xllcorner, self.yllcorner)

    ##
    # @brief Get the X coordinate of the origin of this georef.
    def get_xllcorner(self):
        return self.xllcorner

    ##
    # @brief Get the Y coordinate of the origin of this georef.
    def get_yllcorner(self):
        return self.yllcorner

    ##
    # @brief Get the zone of this georef.
    def get_zone(self):
        return self.zone

    ##
    # @brief Write <something> to an open NetCDF file.
    # @param outfile Handle to open NetCDF file.
    def write_NetCDF(self, outfile):
        outfile.xllcorner = self.xllcorner
        outfile.yllcorner = self.yllcorner
        outfile.zone = self.zone

        outfile.false_easting = self.false_easting
        outfile.false_northing = self.false_northing

        outfile.datum = self.datum
        outfile.projection = self.projection
        outfile.units = self.units

    ##
    # @brief Read data from an open NetCDF file.
    # @param infile Handle to open NetCDF file.
    def read_NetCDF(self, infile):
        self.xllcorner = infile.xllcorner[0]
        self.yllcorner = infile.yllcorner[0]
        self.zone = infile.zone[0]

        # Fix some assertion failures
        if isinstance(self.zone, num.ndarray) and self.zone.shape == ():
            self.zone = self.zone[0]
        if (isinstance(self.xllcorner, num.ndarray) and
                self.xllcorner.shape == ()):
            self.xllcorner = self.xllcorner[0]
        if (isinstance(self.yllcorner, num.ndarray) and
                self.yllcorner.shape == ()):
            self.yllcorner = self.yllcorner[0]

        assert (self.xllcorner.dtype.kind in num.typecodes['Float'] or
                self.xllcorner.dtype.kind in num.typecodes['Integer'])
        assert (self.yllcorner.dtype.kind in num.typecodes['Float'] or
                self.yllcorner.dtype.kind in num.typecodes['Integer'])
        assert (self.zone.dtype.kind in num.typecodes['Integer'])

        try:
            self.false_easting = infile.false_easting[0]
            self.false_northing = infile.false_northing[0]

            self.datum = infile.datum
            self.projection = infile.projection
            self.units = infile.units
        except:
            pass

        if self.false_easting != DEFAULT_FALSE_EASTING:
            print "WARNING: False easting of %f specified." % self.false_easting
            print "Default false easting is %f." % DEFAULT_FALSE_EASTING
            print "ANUGA does not correct for differences in False Eastings."

        if self.false_northing != DEFAULT_FALSE_NORTHING:
            print ("WARNING: False northing of %f specified."
                   % self.false_northing)
            print "Default false northing is %f." % DEFAULT_FALSE_NORTHING
            print "ANUGA does not correct for differences in False Northings."

        if self.datum.upper() != DEFAULT_DATUM.upper():
            print "WARNING: Datum of %s specified." % self.datum
            print "Default Datum is %s." % DEFAULT_DATUM
            print "ANUGA does not correct for differences in datums."

        if self.projection.upper() != DEFAULT_PROJECTION.upper():
            print "WARNING: Projection of %s specified." % self.projection
            print "Default Projection is %s." % DEFAULT_PROJECTION
            print "ANUGA does not correct for differences in Projection."

        if self.units.upper() != DEFAULT_UNITS.upper():
            print "WARNING: Units of %s specified." % self.units
            print "Default units is %s." % DEFAULT_UNITS
            print "ANUGA does not correct for differences in units."

################################################################################
# ASCII files with geo-refs are currently not used
################################################################################

    ##
    # @brief Write georef data to an open text file.
    # @param fd Handle to open text file.
    def write_ASCII(self, fd):
        fd.write(TITLE)
        fd.write(str(self.zone) + "\n")
        fd.write(str(self.xllcorner) + "\n")
        fd.write(str(self.yllcorner) + "\n")

    ##
    # @brief Read georef data from an open text file.
    # @param fd Handle to open text file.
    def read_ASCII(self, fd, read_title=None):
        try:
            if read_title == None:
                read_title = fd.readline()     # remove the title line
            if read_title[0:2].upper() != TITLE[0:2].upper():
                msg = ('File error.  Expecting line: %s.  Got this line: %s'
                       % (TITLE, read_title))
                raise TitleError, msg
            self.zone = int(fd.readline())
            self.xllcorner = float(fd.readline())
            self.yllcorner = float(fd.readline())
        except SyntaxError:
            msg = 'File error.  Got syntax error while parsing geo reference'
            raise ParsingError, msg

        # Fix some assertion failures
        if isinstance(self.zone, num.ndarray) and self.zone.shape == ():
            self.zone = self.zone[0]
        if (isinstance(self.xllcorner, num.ndarray) and
                self.xllcorner.shape == ()):
            self.xllcorner = self.xllcorner[0]
        if (isinstance(self.yllcorner, num.ndarray) and
                self.yllcorner.shape == ()):
            self.yllcorner = self.yllcorner[0]

################################################################################

    def change_points_geo_ref(self, points, points_geo_ref=None):
        """
        Change the geo reference of a list or numeric array of points to
        be this reference.(The reference used for this object)
        If the points do not have a geo ref, assume 'absolute' values
        """

        is_list = False
        if type(points) == types.ListType:
            is_list = True

        points = ensure_numeric(points, num.float)

        if len(points.shape) == 1:
            # One point has been passed
            msg = 'Single point must have two elements'
            assert len(points) == 2, msg
            points = num.reshape(points, (1,2))

        msg = ('Points array must be two dimensional.\n'
               'I got %d dimensions' % len(points.shape))
        assert len(points.shape) == 2, msg

        msg = ('Input must be an N x 2 array or list of (x,y) values. '
               'I got an %d x %d array' % points.shape)
        assert points.shape[1] == 2, msg

        # FIXME (Ole): Could also check if zone, xllcorner, yllcorner
        # are identical in the two geo refs.
        if points_geo_ref is not self:
            # If georeferences are different
            if not points_geo_ref is None:
                # Convert points to absolute coordinates
                points[:,0] += points_geo_ref.xllcorner
                points[:,1] += points_geo_ref.yllcorner

            # Make points relative to primary geo reference
            points[:,0] -= self.xllcorner
            points[:,1] -= self.yllcorner

        # if we got a list, return a list
        if is_list:
            points = points.tolist()

        return points

    ##
    # @brief Is this georef absolute?
    # @return True if ???
    def is_absolute(self):
        """Return True if xllcorner==yllcorner==0
        indicating that points in question are absolute.
        """

        return num.allclose([self.xllcorner, self.yllcorner], 0)

    ##
    # @brief Convert points to absolute points in this georef instance.
    # @param points
    # @return
    # @note
    def get_absolute(self, points):
        """Given a set of points geo referenced to this instance,
        return the points as absolute values.
        """

        #if self.is_absolute:
        #    return points

        # remember if we got a list
        is_list = False
        if type(points) == types.ListType:
            is_list = True

        points = ensure_numeric(points, num.float)
        if len(points.shape) == 1:
            #One point has been passed
            if not len(points) == 2:
                msg = 'Single point must have two elements'
                raise ShapeError, msg
                #points = reshape(points, (1,2))

        if not points.shape[1] == 2:
            msg = ('Input must be an N x 2 array or list of (x,y) values. '
                   'I got an %d x %d array' % points.shape)
            raise ShapeError, msg

        # Add geo ref to points
        #if not self.is_absolute:
        if not self.is_absolute():
            points[:,0] += self.xllcorner
            points[:,1] += self.yllcorner
            #self.is_absolute = True

        # if we got a list, return a list
        if is_list:
            points = points.tolist()

        return points

    ##
    # @brief Convert points to relative measurement.
    # @param points Points to convert to relative measurements.
    # @return A set of points relative to the geo_reference instance.
    def get_relative(self, points):
        """Given a set of points in absolute UTM coordinates,
        make them relative to this geo_reference instance,
        return the points as relative values.

        This is the inverse of get_absolute.
        """

        is_list = False
        if type(points) == types.ListType:
            is_list = True

        points = ensure_numeric(points, num.float)
        if len(points.shape) == 1:
            #One point has been passed
            if not len(points) == 2:
                msg = 'Single point must have two elements'
                raise ShapeError, msg
                #points = reshape(points, (1,2))

        if not points.shape[1] == 2:
            msg = ('Input must be an N x 2 array or list of (x,y) values. '
                   'I got an %d x %d array' % points.shape)
            raise ShapeError, msg

        # Subtract geo ref from points
        if not self.is_absolute():
            points[:,0] -= self.xllcorner
            points[:,1] -= self.yllcorner

        # if we got a list, return a list
        if is_list:
            points = points.tolist()

        return points

    ##
    # @brief ??
    # @param other ??
    def reconcile_zones(self, other):
        if other is None:
            other = Geo_reference()
        if (self.zone == other.zone or
            self.zone == DEFAULT_ZONE and
            other.zone == DEFAULT_ZONE):
            pass
        elif self.zone == DEFAULT_ZONE:
            self.zone = other.zone
        elif other.zone == DEFAULT_ZONE:
            other.zone = self.zone
        else:
            msg = ('Geospatial data must be in the same '
                   'ZONE to allow reconciliation. I got zone %d and %d'
                   % (self.zone, other.zone))
            raise ANUGAError, msg

    #def easting_northing2geo_reffed_point(self, x, y):
    #    return [x-self.xllcorner, y - self.xllcorner]

    #def easting_northing2geo_reffed_points(self, x, y):
    #    return [x-self.xllcorner, y - self.xllcorner]

    ##
    # @brief Get origin of this geo_reference.
    # @return (zone, xllcorner, yllcorner).
    def get_origin(self):
        return (self.zone, self.xllcorner, self.yllcorner)

    ##
    # @brief Get a string representation of this geo_reference instance.
    def __repr__(self):
        return ('(zone=%i easting=%f, northing=%f)'
                % (self.zone, self.xllcorner, self.yllcorner))

    ##
    # @brief Compare two geo_reference instances.
    # @param self This geo_reference instance.
    # @param other Another geo_reference instance to compare against.
    # @return 0 if instances have the same attributes, else 1.
    # @note Attributes are: zone, xllcorner, yllcorner.
    def __cmp__(self, other):
        # FIXME (DSG) add a tolerence
        if other is None:
            return 1
        cmp = 0
        if not (self.xllcorner == self.xllcorner):
            cmp = 1
        if not (self.yllcorner == self.yllcorner):
            cmp = 1
        if not (self.zone == self.zone):
            cmp = 1
        return cmp


##
# @brief Write a geo_reference to a NetCDF file (usually SWW).
# @param origin A georef instance or parameters to create a georef instance.
# @param outfile Path to file to write.
# @return A normalized geo_reference.
def write_NetCDF_georeference(origin, outfile):
    """Write georeference info to a netcdf file, usually sww.

    The origin can be a georef instance or parameters for a geo_ref instance

    outfile is the name of the file to be written to.
    """

    geo_ref = ensure_geo_reference(origin)
    geo_ref.write_NetCDF(outfile)
    return geo_ref


##
# @brief Convert an object to a georeference instance.
# @param origin A georef instance or (zone, xllcorner, yllcorner)
# @return A georef object, or None if 'origin' was None.
def ensure_geo_reference(origin):
    """
    Given a list/tuple of zone, xllcorner and yllcorner of a geo-ref object,
    return a geo ref object.

    If the origin is None, return None, so calling this function doesn't
    effect code logic
    """

    if isinstance(origin, Geo_reference):
        geo_ref = origin
    elif origin is None:
        geo_ref = None
    else:
        geo_ref = apply(Geo_reference, origin)

    return geo_ref


#-----------------------------------------------------------------------

if __name__ == "__main__":
    pass