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

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

Added tests for get_relative().

File size: 16.6 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        #self.is_absolute = num.allclose([self.xllcorner, self.yllcorner], 0)
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    ##
106    # @brief Get the X coordinate of the origin of this georef.
107    def get_xllcorner(self):
108        return self.xllcorner
109
110    ##
111    # @brief Get the Y coordinate of the origin of this georef.
112    def get_yllcorner(self):
113        return self.yllcorner
114
115    ##
116    # @brief Get the zone of this georef.
117    def get_zone(self):
118        return self.zone
119
120    ##
121    # @brief Write <something> to an open NetCDF file.
122    # @param outfile Handle to open NetCDF file.
123    def write_NetCDF(self, outfile):
124        outfile.xllcorner = self.xllcorner
125        outfile.yllcorner = self.yllcorner
126        outfile.zone = self.zone
127
128        outfile.false_easting = self.false_easting
129        outfile.false_northing = self.false_northing
130
131        outfile.datum = self.datum
132        outfile.projection = self.projection
133        outfile.units = self.units
134
135    ##
136    # @brief Read data from an open NetCDF file.
137    # @param infile Handle to open NetCDF file.
138    def read_NetCDF(self, infile):
139        self.xllcorner = infile.xllcorner[0]
140        self.yllcorner = infile.yllcorner[0]
141        self.zone = infile.zone[0]
142
143        # Fix some assertion failures
144        if isinstance(self.zone, num.ndarray) and self.zone.shape == ():
145            self.zone = self.zone[0]
146        if (isinstance(self.xllcorner, num.ndarray) and
147                self.xllcorner.shape == ()):
148            self.xllcorner = self.xllcorner[0]
149        if (isinstance(self.yllcorner, num.ndarray) and
150                self.yllcorner.shape == ()):
151            self.yllcorner = self.yllcorner[0]
152
153        assert (self.xllcorner.dtype.kind in num.typecodes['Float'] or
154                self.xllcorner.dtype.kind in num.typecodes['Integer'])
155        assert (self.yllcorner.dtype.kind in num.typecodes['Float'] or
156                self.yllcorner.dtype.kind in num.typecodes['Integer'])
157        assert (self.zone.dtype.kind in num.typecodes['Integer'])
158
159        try:
160            self.false_easting = infile.false_easting[0]
161            self.false_northing = infile.false_northing[0]
162
163            self.datum = infile.datum
164            self.projection = infile.projection
165            self.units = infile.units
166        except:
167            pass
168
169        if self.false_easting != DEFAULT_FALSE_EASTING:
170            print "WARNING: False easting of %f specified." % self.false_easting
171            print "Default false easting is %f." % DEFAULT_FALSE_EASTING
172            print "ANUGA does not correct for differences in False Eastings."
173
174        if self.false_northing != DEFAULT_FALSE_NORTHING:
175            print ("WARNING: False northing of %f specified."
176                   % self.false_northing)
177            print "Default false northing is %f." % DEFAULT_FALSE_NORTHING
178            print "ANUGA does not correct for differences in False Northings."
179
180        if self.datum.upper() != DEFAULT_DATUM.upper():
181            print "WARNING: Datum of %s specified." % self.datum
182            print "Default Datum is %s." % DEFAULT_DATUM
183            print "ANUGA does not correct for differences in datums."
184
185        if self.projection.upper() != DEFAULT_PROJECTION.upper():
186            print "WARNING: Projection of %s specified." % self.projection
187            print "Default Projection is %s." % DEFAULT_PROJECTION
188            print "ANUGA does not correct for differences in Projection."
189
190        if self.units.upper() != DEFAULT_UNITS.upper():
191            print "WARNING: Units of %s specified." % self.units
192            print "Default units is %s." % DEFAULT_UNITS
193            print "ANUGA does not correct for differences in units."
194
195################################################################################
196# ASCII files with geo-refs are currently not used
197################################################################################
198
199    ##
200    # @brief Write georef data to an open text file.
201    # @param fd Handle to open text file.
202    def write_ASCII(self, fd):
203        fd.write(TITLE)
204        fd.write(str(self.zone) + "\n")
205        fd.write(str(self.xllcorner) + "\n")
206        fd.write(str(self.yllcorner) + "\n")
207
208    ##
209    # @brief Read georef data from an open text file.
210    # @param fd Handle to open text file.
211    def read_ASCII(self, fd, read_title=None):
212        try:
213            if read_title == None:
214                read_title = fd.readline()     # remove the title line
215            if read_title[0:2].upper() != TITLE[0:2].upper():
216                msg = ('File error.  Expecting line: %s.  Got this line: %s'
217                       % (TITLE, read_title))
218                raise TitleError, msg
219            self.zone = int(fd.readline())
220            self.xllcorner = float(fd.readline())
221            self.yllcorner = float(fd.readline())
222        except SyntaxError:
223            msg = 'File error.  Got syntax error while parsing geo reference'
224            raise ParsingError, msg
225
226        # Fix some assertion failures
227        if isinstance(self.zone, num.ndarray) and self.zone.shape == ():
228            self.zone = self.zone[0]
229        if (isinstance(self.xllcorner, num.ndarray) and
230                self.xllcorner.shape == ()):
231            self.xllcorner = self.xllcorner[0]
232        if (isinstance(self.yllcorner, num.ndarray) and
233                self.yllcorner.shape == ()):
234            self.yllcorner = self.yllcorner[0]
235
236################################################################################
237
238    ##
239    # @brief Change points to be absolute wrt new georef 'points_geo_ref'.
240    # @param points The points to change.
241    # @param points_geo_ref The new georef to make points absolute wrt.
242    # @return The changed points.
243    # @note If 'points' is a list then a changed list is returned.
244    def change_points_geo_ref(self, points, points_geo_ref=None):
245        """Change the geo reference of a list or numeric array of points to
246        be this reference.(The reference used for this object)
247        If the points do not have a geo ref, assume 'absolute' input values
248        """
249
250        # remember if we got a list
251        is_list = isinstance(points, list)
252
253        points = ensure_numeric(copy.copy(points), num.float)
254
255        # sanity checks
256        if len(points.shape) == 1:
257            # One point has been passed
258            msg = 'Single point must have two elements'
259            assert len(points) == 2, msg
260            points = num.reshape(points, (1,2))
261
262        msg = ('Points array must be two dimensional.\n'
263               'I got %d dimensions' % len(points.shape))
264        assert len(points.shape) == 2, msg
265
266        msg = ('Input must be an N x 2 array or list of (x,y) values. '
267               'I got an %d x %d array' % points.shape)
268        assert points.shape[1] == 2, msg
269
270        # FIXME (Ole): Could also check if zone, xllcorner, yllcorner
271        # are identical in the two geo refs.
272        if points_geo_ref is not self:
273            # If georeferences are different
274            if not points_geo_ref is None:
275                # Convert points to absolute coordinates
276                points[:,0] += points_geo_ref.xllcorner
277                points[:,1] += points_geo_ref.yllcorner
278
279            # Make points relative to primary geo reference
280            points[:,0] -= self.xllcorner
281            points[:,1] -= self.yllcorner
282
283        # if we got a list, return a list
284        if is_list:
285            points = points.tolist()
286
287        return points
288
289    ##
290    # @brief Is this georef absolute?
291    # @return True if ???
292    def is_absolute(self):
293        """Return True if xllcorner==yllcorner==0"""
294
295        return num.allclose([self.xllcorner, self.yllcorner], 0)
296
297    ##
298    # @brief Convert points to absolute points in this georef instance.
299    # @param points
300    # @return
301    # @note This method changes the input points!
302    def get_absolute(self, points):
303        """Given a set of points geo referenced to this instance
304
305        return the points as absolute values.
306        """
307
308        # remember if we got a list
309        is_list = isinstance(points, list)
310
311        # convert to numeric array, force a copy
312        points = ensure_numeric(copy.copy(points), num.float)
313
314        # sanity checks
315        if len(points.shape) == 1:
316            #One point has been passed
317            if not len(points) == 2:
318                msg = 'Single point must have two elements'
319                raise ShapeError, msg
320
321        if not points.shape[1] == 2:
322            msg = ('Input must be an N x 2 array or list of (x,y) values. '
323                   'I got an %d x %d array' % points.shape)
324            raise ShapeError, msg
325
326        # Add geo ref to points
327        if not self.is_absolute():
328            points[:,0] += self.xllcorner
329            points[:,1] += self.yllcorner
330
331        # if we got a list, return a list
332        if is_list:
333            points = points.tolist()
334
335        return points
336
337    ##
338    # @brief Convert points to relative measurement.
339    # @param points Points to convert to relative measurements.
340    # @return A set of points relative to the geo_reference instance.
341    def get_relative(self, points):
342        """Given a set of points in absolute UTM coordinates,
343        make them relative to this geo_reference instance,
344        return the points as relative values.
345
346        This is the inverse of get_absolute.
347        """
348
349        # remember if we got a list
350        is_list = isinstance(points, list)
351
352        # convert to a numeric array, force a copy
353        points = ensure_numeric(copy.copy(points), num.float)
354
355        # sanity checks
356        if len(points.shape) == 1:
357            #One point has been passed
358            if not len(points) == 2:
359                msg = 'Single point must have two elements'
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[:,0] -= self.xllcorner
370            points[:,1] -= self.yllcorner
371
372        # if we got a list, return a list
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.