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

Last change on this file since 7035 was 7035, checked in by rwilson, 15 years ago

Back-merge from Numeric trunk.

File size: 17.4 KB
Line 
1"""
2
3"""
4
5
6#FIXME: Ensure that all attributes of a georef are treated everywhere
7#and unit test
8
9import types, sys
10import copy
11
12from anuga.utilities.numerical_tools import ensure_numeric
13from anuga.utilities.anuga_exceptions import ANUGAError, TitleError, \
14                                             ParsingError, ShapeError
15from anuga.config import netcdf_float, netcdf_int, netcdf_float32
16
17import numpy as num
18
19
20DEFAULT_ZONE = -1
21TITLE = '#geo reference' + "\n" # this title is referred to in the test format
22
23DEFAULT_PROJECTION = 'UTM'
24DEFAULT_DATUM = 'wgs84'
25DEFAULT_UNITS = 'm'
26DEFAULT_FALSE_EASTING = 500000
27DEFAULT_FALSE_NORTHING = 10000000    # Default for southern hemisphere
28
29
30##
31# @brief A class for ...
32class Geo_reference:
33    """
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
45    """
46
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.
60    def __init__(self,
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,
69                 NetCDFObject=None,
70                 ASCIIFile=None,
71                 read_title=None):
72        """
73        input:
74        NetCDFObject - a handle to the netCDF file to be written to
75        ASCIIFile - a handle to the text file
76        read_title - the title of the georeference text, if it was read in.
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.
80
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.
85        """
86
87        if zone is None:
88            zone = DEFAULT_ZONE
89        self.false_easting = false_easting
90        self.false_northing = false_northing
91        self.datum = datum
92        self.projection = projection
93        self.zone = zone
94        self.units = units
95        self.xllcorner = xllcorner
96        self.yllcorner = yllcorner       
97           
98        if NetCDFObject is not None:
99            self.read_NetCDF(NetCDFObject)
100
101        if ASCIIFile is not None:
102            self.read_ASCII(ASCIIFile, read_title=read_title)
103
104           
105        # Set flag for absolute points (used by get_absolute)   
106        self.absolute = num.allclose([self.xllcorner, self.yllcorner], 0)
107           
108
109    def get_xllcorner(self):
110        return self.xllcorner
111
112    ##
113    # @brief Get the Y coordinate of the origin of this georef.
114    def get_yllcorner(self):
115        return self.yllcorner
116
117    ##
118    # @brief Get the zone of this georef.
119    def get_zone(self):
120        return self.zone
121
122    ##
123    # @brief Write <something> to an open NetCDF file.
124    # @param outfile Handle to open NetCDF file.
125    def write_NetCDF(self, outfile):
126        outfile.xllcorner = self.xllcorner
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
133        outfile.datum = self.datum
134        outfile.projection = self.projection
135        outfile.units = self.units
136
137    ##
138    # @brief Read data from an open NetCDF file.
139    # @param infile Handle to open NetCDF file.
140    def read_NetCDF(self, infile):
141        self.xllcorner = infile.xllcorner[0]
142        self.yllcorner = infile.yllcorner[0]
143        self.zone = infile.zone[0]
144
145        # Fix some assertion failures
146        if isinstance(self.zone, num.ndarray) and self.zone.shape == ():
147            self.zone = self.zone[0]
148        if (isinstance(self.xllcorner, num.ndarray) and
149                self.xllcorner.shape == ()):
150            self.xllcorner = self.xllcorner[0]
151        if (isinstance(self.yllcorner, num.ndarray) and
152                self.yllcorner.shape == ()):
153            self.yllcorner = self.yllcorner[0]
154
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'])
160
161        try:
162            self.false_easting = infile.false_easting[0]
163            self.false_northing = infile.false_northing[0]
164
165            self.datum = infile.datum
166            self.projection = infile.projection
167            self.units = infile.units
168        except:
169            pass
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
174            print "ANUGA does not correct for differences in False Eastings."
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
180            print "ANUGA does not correct for differences in False Northings."
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
185            print "ANUGA does not correct for differences in datums."
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
190            print "ANUGA does not correct for differences in Projection."
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
195            print "ANUGA does not correct for differences in units."
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.
204    def write_ASCII(self, fd):
205        fd.write(TITLE)
206        fd.write(str(self.zone) + "\n")
207        fd.write(str(self.xllcorner) + "\n")
208        fd.write(str(self.yllcorner) + "\n")
209
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):
214        try:
215            if read_title == None:
216                read_title = fd.readline()     # remove the title line
217            if read_title[0:2].upper() != TITLE[0:2].upper():
218                msg = ('File error.  Expecting line: %s.  Got this line: %s'
219                       % (TITLE, read_title))
220                raise TitleError, msg
221            self.zone = int(fd.readline())
222            self.xllcorner = float(fd.readline())
223            self.yllcorner = float(fd.readline())
224        except SyntaxError:
225            msg = 'File error.  Got syntax error while parsing geo reference'
226            raise ParsingError, msg
227
228        # Fix some assertion failures
229        if isinstance(self.zone, num.ndarray) and self.zone.shape == ():
230            self.zone = self.zone[0]
231        if (isinstance(self.xllcorner, num.ndarray) and
232                self.xllcorner.shape == ()):
233            self.xllcorner = self.xllcorner[0]
234        if (isinstance(self.yllcorner, num.ndarray) and
235                self.yllcorner.shape == ()):
236            self.yllcorner = self.yllcorner[0]
237
238        assert (type(self.xllcorner) == types.FloatType)
239        assert (type(self.yllcorner) == types.FloatType)
240        assert (type(self.zone) == types.IntType)
241
242################################################################################
243
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.
250    def change_points_geo_ref(self, points, points_geo_ref=None):
251        """Change the geo reference of a list or Numeric array of points to
252        be this reference.(The reference used for this object)
253        If the points do not have a geo ref, assume 'absolute' values
254        """
255        import copy
256       
257        # remember if we got a list
258        is_list = isinstance(points, list)
259
260        points = ensure_numeric(points, num.float)
261
262        # sanity checks
263        if len(points.shape) == 1:
264            #One point has been passed
265            msg = 'Single point must have two elements'
266            assert len(points) == 2, msg
267            points = num.reshape(points, (1,2))
268
269        msg = 'Points array must be two dimensional.\n'
270        msg += 'I got %d dimensions' %len(points.shape)
271        assert len(points.shape) == 2, msg
272
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               
276
277        # FIXME (Ole): Could also check if zone, xllcorner, yllcorner
278        # are identical in the two geo refs.   
279        if points_geo_ref is not self:
280            # If georeferences are different
281            points = copy.copy(points) # Don't destroy input                   
282            if not points_geo_ref is None:
283                # Convert points to absolute coordinates
284                points[:,0] += points_geo_ref.xllcorner
285                points[:,1] += points_geo_ref.yllcorner
286       
287            # Make points relative to primary geo reference
288            points[:,0] -= self.xllcorner
289            points[:,1] -= self.yllcorner
290
291        if is_list:
292            points = points.tolist()
293           
294        return points
295
296    def is_absolute(self):
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).
307
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
315
316    def get_absolute(self, points):
317        """Given a set of points geo referenced to this instance,
318        return the points as absolute values.
319        """
320
321        # remember if we got a list
322        is_list = isinstance(points, list)
323
324        points = ensure_numeric(points, num.float)
325        if len(points.shape) == 1:
326            # One point has been passed
327            msg = 'Single point must have two elements'
328            if not len(points) == 2:
329                raise ShapeError, msg   
330
331
332        msg = 'Input must be an N x 2 array or list of (x,y) values. '
333        msg += 'I got an %d x %d array' %points.shape   
334        if not points.shape[1] == 2:
335            raise ShapeError, msg   
336           
337       
338        # Add geo ref to points
339        if not self.is_absolute():
340            points = copy.copy(points) # Don't destroy input                   
341            points[:,0] += self.xllcorner
342            points[:,1] += self.yllcorner
343
344       
345        if is_list:
346            points = points.tolist()
347             
348        return points
349
350    ##
351    # @brief Convert points to relative measurement.
352    # @param points Points to convert to relative measurements.
353    # @return A set of points relative to the geo_reference instance.
354    def get_relative(self, points):
355        """Given a set of points in absolute UTM coordinates,
356        make them relative to this geo_reference instance,
357        return the points as relative values.
358
359        This is the inverse of get_absolute.
360        """
361
362        # remember if we got a list
363        is_list = isinstance(points, list)
364
365        points = ensure_numeric(points, num.float)
366        if len(points.shape) == 1:
367            #One point has been passed
368            msg = 'Single point must have two elements'
369            if not len(points) == 2:
370                raise ShapeError, msg   
371
372        if not points.shape[1] == 2:
373            msg = ('Input must be an N x 2 array or list of (x,y) values. '
374                   'I got an %d x %d array' % points.shape)
375            raise ShapeError, msg   
376
377        # Subtract geo ref from points
378        if not self.is_absolute():
379            points = copy.copy(points) # Don't destroy input                           
380            points[:,0] -= self.xllcorner
381            points[:,1] -= self.yllcorner
382
383        if is_list:
384            points = points.tolist()
385             
386        return points
387
388    ##
389    # @brief ??
390    # @param other ??
391    def reconcile_zones(self, other):
392        if other is None:
393            other = Geo_reference()
394        if (self.zone == other.zone or
395            self.zone == DEFAULT_ZONE and
396            other.zone == DEFAULT_ZONE):
397            pass
398        elif self.zone == DEFAULT_ZONE:
399            self.zone = other.zone
400        elif other.zone == DEFAULT_ZONE:
401            other.zone = self.zone
402        else:
403            msg = ('Geospatial data must be in the same '
404                   'ZONE to allow reconciliation. I got zone %d and %d'
405                   % (self.zone, other.zone))
406            raise ANUGAError, msg
407
408    #def easting_northing2geo_reffed_point(self, x, y):
409    #    return [x-self.xllcorner, y - self.xllcorner]
410
411    #def easting_northing2geo_reffed_points(self, x, y):
412    #    return [x-self.xllcorner, y - self.xllcorner]
413
414    ##
415    # @brief Get origin of this geo_reference.
416    # @return (zone, xllcorner, yllcorner).
417    def get_origin(self):
418        return (self.zone, self.xllcorner, self.yllcorner)
419
420    ##
421    # @brief Get a string representation of this geo_reference instance.
422    def __repr__(self):
423        return ('(zone=%i easting=%f, northing=%f)'
424                % (self.zone, self.xllcorner, self.yllcorner))
425
426    ##
427    # @brief Compare two geo_reference instances.
428    # @param self This geo_reference instance.
429    # @param other Another geo_reference instance to compare against.
430    # @return 0 if instances have the same attributes, else 1.
431    # @note Attributes are: zone, xllcorner, yllcorner.
432    def __cmp__(self, other):
433        # FIXME (DSG) add a tolerence
434        if other is None:
435            return 1
436        cmp = 0
437        if not (self.xllcorner == self.xllcorner):
438            cmp = 1
439        if not (self.yllcorner == self.yllcorner):
440            cmp = 1
441        if not (self.zone == self.zone):
442            cmp = 1
443        return cmp
444
445
446##
447# @brief Write a geo_reference to a NetCDF file (usually SWW).
448# @param origin A georef instance or parameters to create a georef instance.
449# @param outfile Path to file to write.
450# @return A normalized geo_reference.
451def write_NetCDF_georeference(origin, outfile):
452    """Write georeference info to a netcdf file, usually sww.
453
454    The origin can be a georef instance or parameters for a geo_ref instance
455
456    outfile is the name of the file to be written to.
457    """
458
459    geo_ref = ensure_geo_reference(origin)
460    geo_ref.write_NetCDF(outfile)
461    return geo_ref
462
463
464##
465# @brief Convert an object to a georeference instance.
466# @param origin A georef instance or (zone, xllcorner, yllcorner)
467# @return A georef object, or None if 'origin' was None.
468def ensure_geo_reference(origin):
469    """
470    Given a list/tuple of zone, xllcorner and yllcorner of a geo-ref object,
471    return a geo ref object.
472
473    If the origin is None, return None, so calling this function doesn't
474    effect code logic
475    """
476
477    if isinstance(origin, Geo_reference):
478        geo_ref = origin
479    elif origin is None:
480        geo_ref = None
481    else:
482        geo_ref = apply(Geo_reference, origin)
483
484    return geo_ref
485
486
487#-----------------------------------------------------------------------
488
489if __name__ == "__main__":
490    pass
Note: See TracBrowser for help on using the repository browser.