source: trunk/anuga_core/source/anuga/coordinate_transforms/geo_reference.py @ 7967

Last change on this file since 7967 was 7779, checked in by hudson, 14 years ago

Fixed wrong exception module import paths.

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