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

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

Back-merge from Numeric trunk.
File size: 16.7 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 = int(false_easting)
90        self.false_northing = int(false_northing)
91        self.datum = datum
92        self.projection = projection
93        self.zone = int(zone)
94        self.units = units
95        self.xllcorner = float(xllcorner)
96        self.yllcorner = float(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        # Set flag for absolute points (used by get_absolute)   
105        self.absolute = num.allclose([self.xllcorner, self.yllcorner], 0)
106           
107
108    def get_xllcorner(self):
109        return self.xllcorner
110
111    ##
112    # @brief Get the Y coordinate of the origin of this georef.
113    def get_yllcorner(self):
114        return self.yllcorner
115
116    ##
117    # @brief Get the zone of this georef.
118    def get_zone(self):
119        return self.zone
120
121    ##
122    # @brief Write <something> to an open NetCDF file.
123    # @param outfile Handle to open NetCDF file.
124    def write_NetCDF(self, outfile):
125        outfile.xllcorner = self.xllcorner
126        outfile.yllcorner = self.yllcorner
127        outfile.zone = self.zone
128
129        outfile.false_easting = self.false_easting
130        outfile.false_northing = self.false_northing
131
132        outfile.datum = self.datum
133        outfile.projection = self.projection
134        outfile.units = self.units
135
136    ##
137    # @brief Read data from an open NetCDF file.
138    # @param infile Handle to open NetCDF file.
139    def read_NetCDF(self, infile):
140        self.xllcorner = float(infile.xllcorner[0])
141        self.yllcorner = float(infile.yllcorner[0])
142        self.zone = int(infile.zone[0])
143
144        try:
145            self.false_easting = int(infile.false_easting[0])
146            self.false_northing = int(infile.false_northing[0])
147
148            self.datum = infile.datum
149            self.projection = infile.projection
150            self.units = infile.units
151        except:
152            pass
153
154        if self.false_easting != DEFAULT_FALSE_EASTING:
155            print "WARNING: False easting of %f specified." % self.false_easting
156            print "Default false easting is %f." % DEFAULT_FALSE_EASTING
157            print "ANUGA does not correct for differences in False Eastings."
158
159        if self.false_northing != DEFAULT_FALSE_NORTHING:
160            print ("WARNING: False northing of %f specified."
161                   % self.false_northing)
162            print "Default false northing is %f." % DEFAULT_FALSE_NORTHING
163            print "ANUGA does not correct for differences in False Northings."
164
165        if self.datum.upper() != DEFAULT_DATUM.upper():
166            print "WARNING: Datum of %s specified." % self.datum
167            print "Default Datum is %s." % DEFAULT_DATUM
168            print "ANUGA does not correct for differences in datums."
169
170        if self.projection.upper() != DEFAULT_PROJECTION.upper():
171            print "WARNING: Projection of %s specified." % self.projection
172            print "Default Projection is %s." % DEFAULT_PROJECTION
173            print "ANUGA does not correct for differences in Projection."
174
175        if self.units.upper() != DEFAULT_UNITS.upper():
176            print "WARNING: Units of %s specified." % self.units
177            print "Default units is %s." % DEFAULT_UNITS
178            print "ANUGA does not correct for differences in units."
179
180################################################################################
181# ASCII files with geo-refs are currently not used
182################################################################################
183
184    ##
185    # @brief Write georef data to an open text file.
186    # @param fd Handle to open text file.
187    def write_ASCII(self, fd):
188        fd.write(TITLE)
189        fd.write(str(self.zone) + "\n")
190        fd.write(str(self.xllcorner) + "\n")
191        fd.write(str(self.yllcorner) + "\n")
192
193    ##
194    # @brief Read georef data from an open text file.
195    # @param fd Handle to open text file.
196    def read_ASCII(self, fd, read_title=None):
197        try:
198            if read_title == None:
199                read_title = fd.readline()     # remove the title line
200            if read_title[0:2].upper() != TITLE[0:2].upper():
201                msg = ('File error.  Expecting line: %s.  Got this line: %s'
202                       % (TITLE, read_title))
203                raise TitleError, msg
204            self.zone = int(fd.readline())
205            self.xllcorner = float(fd.readline())
206            self.yllcorner = float(fd.readline())
207        except SyntaxError:
208            msg = 'File error.  Got syntax error while parsing geo reference'
209            raise ParsingError, msg
210
211        # Fix some assertion failures
212        if isinstance(self.zone, num.ndarray) and self.zone.shape == ():
213            self.zone = self.zone[0]
214        if (isinstance(self.xllcorner, num.ndarray) and
215                self.xllcorner.shape == ()):
216            self.xllcorner = self.xllcorner[0]
217        if (isinstance(self.yllcorner, num.ndarray) and
218                self.yllcorner.shape == ()):
219            self.yllcorner = self.yllcorner[0]
220
221        assert (type(self.xllcorner) == types.FloatType)
222        assert (type(self.yllcorner) == types.FloatType)
223        assert (type(self.zone) == types.IntType)
224
225################################################################################
226
227    ##
228    # @brief Change points to be absolute wrt new georef 'points_geo_ref'.
229    # @param points The points to change.
230    # @param points_geo_ref The new georef to make points absolute wrt.
231    # @return The changed points.
232    # @note If 'points' is a list then a changed list is returned.
233    def change_points_geo_ref(self, points, points_geo_ref=None):
234        """Change the geo reference of a list or numeric array of points to
235        be this reference.(The reference used for this object)
236        If the points do not have a geo ref, assume 'absolute' values
237        """
238        import copy
239       
240        # remember if we got a list
241        is_list = isinstance(points, list)
242
243        points = ensure_numeric(points, num.float)
244
245        # sanity checks
246        if len(points.shape) == 1:
247            #One point has been passed
248            msg = 'Single point must have two elements'
249            assert len(points) == 2, msg
250            points = num.reshape(points, (1,2))
251
252        msg = 'Points array must be two dimensional.\n'
253        msg += 'I got %d dimensions' %len(points.shape)
254        assert len(points.shape) == 2, msg
255
256        msg = 'Input must be an N x 2 array or list of (x,y) values. '
257        msg += 'I got an %d x %d array' %points.shape   
258        assert points.shape[1] == 2, msg               
259
260        # FIXME (Ole): Could also check if zone, xllcorner, yllcorner
261        # are identical in the two geo refs.   
262        if points_geo_ref is not self:
263            # If georeferences are different
264            points = copy.copy(points) # Don't destroy input                   
265            if not points_geo_ref is None:
266                # Convert points to absolute coordinates
267                points[:,0] += points_geo_ref.xllcorner
268                points[:,1] += points_geo_ref.yllcorner
269       
270            # Make points relative to primary geo reference
271            points[:,0] -= self.xllcorner
272            points[:,1] -= self.yllcorner
273
274        if is_list:
275            points = points.tolist()
276           
277        return points
278
279    def is_absolute(self):
280        """Return True if xllcorner==yllcorner==0 indicating that points
281        in question are absolute.
282        """
283       
284        # FIXME(Ole): It is unfortunate that decision about whether points
285        # are absolute or not lies with the georeference object. Ross pointed this out.
286        # Moreover, this little function is responsible for a large fraction of the time
287        # using in data fitting (something in like 40 - 50%.
288        # This was due to the repeated calls to allclose.
289        # With the flag method fitting is much faster (18 Mar 2009).
290
291        # FIXME(Ole): HACK to be able to reuse data already cached (18 Mar 2009).
292        # Remove at some point
293        if not hasattr(self, 'absolute'):
294            self.absolute = num.allclose([self.xllcorner, self.yllcorner], 0)
295           
296        # Return absolute flag   
297        return self.absolute
298
299    def get_absolute(self, points):
300        """Given a set of points geo referenced to this instance,
301        return the points as absolute values.
302        """
303
304        # remember if we got a list
305        is_list = isinstance(points, list)
306
307        points = ensure_numeric(points, num.float)
308        if len(points.shape) == 1:
309            # One point has been passed
310            msg = 'Single point must have two elements'
311            if not len(points) == 2:
312                raise ShapeError, msg   
313
314
315        msg = 'Input must be an N x 2 array or list of (x,y) values. '
316        msg += 'I got an %d x %d array' %points.shape   
317        if not points.shape[1] == 2:
318            raise ShapeError, msg   
319           
320       
321        # Add geo ref to points
322        if not self.is_absolute():
323            points = copy.copy(points) # Don't destroy input                   
324            points[:,0] += self.xllcorner
325            points[:,1] += self.yllcorner
326
327       
328        if is_list:
329            points = points.tolist()
330             
331        return points
332
333    ##
334    # @brief Convert points to relative measurement.
335    # @param points Points to convert to relative measurements.
336    # @return A set of points relative to the geo_reference instance.
337    def get_relative(self, points):
338        """Given a set of points in absolute UTM coordinates,
339        make them relative to this geo_reference instance,
340        return the points as relative values.
341
342        This is the inverse of get_absolute.
343        """
344
345        # remember if we got a list
346        is_list = isinstance(points, list)
347
348        points = ensure_numeric(points, num.float)
349        if len(points.shape) == 1:
350            #One point has been passed
351            msg = 'Single point must have two elements'
352            if not len(points) == 2:
353                raise ShapeError, msg   
354
355        if not points.shape[1] == 2:
356            msg = ('Input must be an N x 2 array or list of (x,y) values. '
357                   'I got an %d x %d array' % points.shape)
358            raise ShapeError, msg   
359
360        # Subtract geo ref from points
361        if not self.is_absolute():
362            points = copy.copy(points) # Don't destroy input                           
363            points[:,0] -= self.xllcorner
364            points[:,1] -= self.yllcorner
365
366        if is_list:
367            points = points.tolist()
368             
369        return points
370
371    ##
372    # @brief ??
373    # @param other ??
374    def reconcile_zones(self, other):
375        if other is None:
376            other = Geo_reference()
377        if (self.zone == other.zone or
378            self.zone == DEFAULT_ZONE and
379            other.zone == DEFAULT_ZONE):
380            pass
381        elif self.zone == DEFAULT_ZONE:
382            self.zone = other.zone
383        elif other.zone == DEFAULT_ZONE:
384            other.zone = self.zone
385        else:
386            msg = ('Geospatial data must be in the same '
387                   'ZONE to allow reconciliation. I got zone %d and %d'
388                   % (self.zone, other.zone))
389            raise ANUGAError, msg
390
391    #def easting_northing2geo_reffed_point(self, x, y):
392    #    return [x-self.xllcorner, y - self.xllcorner]
393
394    #def easting_northing2geo_reffed_points(self, x, y):
395    #    return [x-self.xllcorner, y - self.xllcorner]
396
397    ##
398    # @brief Get origin of this geo_reference.
399    # @return (zone, xllcorner, yllcorner).
400    def get_origin(self):
401        return (self.zone, self.xllcorner, self.yllcorner)
402
403    ##
404    # @brief Get a string representation of this geo_reference instance.
405    def __repr__(self):
406        return ('(zone=%i easting=%f, northing=%f)'
407                % (self.zone, self.xllcorner, self.yllcorner))
408
409    ##
410    # @brief Compare two geo_reference instances.
411    # @param self This geo_reference instance.
412    # @param other Another geo_reference instance to compare against.
413    # @return 0 if instances have the same attributes, else 1.
414    # @note Attributes are: zone, xllcorner, yllcorner.
415    def __cmp__(self, other):
416        # FIXME (DSG) add a tolerence
417        if other is None:
418            return 1
419        cmp = 0
420        if not (self.xllcorner == self.xllcorner):
421            cmp = 1
422        if not (self.yllcorner == self.yllcorner):
423            cmp = 1
424        if not (self.zone == self.zone):
425            cmp = 1
426        return cmp
427
428
429##
430# @brief Write a geo_reference to a NetCDF file (usually SWW).
431# @param origin A georef instance or parameters to create a georef instance.
432# @param outfile Path to file to write.
433# @return A normalized geo_reference.
434def write_NetCDF_georeference(origin, outfile):
435    """Write georeference info to a netcdf file, usually sww.
436
437    The origin can be a georef instance or parameters for a geo_ref instance
438
439    outfile is the name of the file to be written to.
440    """
441
442    geo_ref = ensure_geo_reference(origin)
443    geo_ref.write_NetCDF(outfile)
444    return geo_ref
445
446
447##
448# @brief Convert an object to a georeference instance.
449# @param origin A georef instance or (zone, xllcorner, yllcorner)
450# @return A georef object, or None if 'origin' was None.
451def ensure_geo_reference(origin):
452    """
453    Given a list/tuple of zone, xllcorner and yllcorner of a geo-ref object,
454    return a geo ref object.
455
456    If the origin is None, return None, so calling this function doesn't
457    effect code logic
458    """
459
460    if isinstance(origin, Geo_reference):
461        geo_ref = origin
462    elif origin is None:
463        geo_ref = None
464    else:
465        geo_ref = apply(Geo_reference, origin)
466
467    return geo_ref
468
469
470#-----------------------------------------------------------------------
471
472if __name__ == "__main__":
473    pass
Note: See TracBrowser for help on using the repository browser.