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

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

Ongoing conversion changes.

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
10from anuga.utilities.numerical_tools import ensure_numeric
11from anuga.utilities.anuga_exceptions import ANUGAError, TitleError, \
12                                             ParsingError, ShapeError
13from anuga.config import netcdf_float, netcdf_int, netcdf_float32
14
15import numpy as num
16
17
18DEFAULT_ZONE = -1
19TITLE = '#geo reference' + "\n" # this title is referred to in the test format
20
21DEFAULT_PROJECTION = 'UTM'
22DEFAULT_DATUM = 'wgs84'
23DEFAULT_UNITS = 'm'
24DEFAULT_FALSE_EASTING = 500000
25DEFAULT_FALSE_NORTHING = 10000000    # Default for southern hemisphere
26
27
28##
29# @brief A class for ...
30class Geo_reference:
31    """
32    Attributes of the Geo_reference class:
33        .zone           The UTM zone (default is -1)
34        .false_easting  ??
35        .false_northing ??
36        .datum          The Datum used (default is wgs84)
37        .projection     The projection used (default is 'UTM')
38        .units          The units of measure used (default metres)
39        .xllcorner      The X coord of origin (default is 0.0 wrt UTM grid)
40        .yllcorner      The y coord of origin (default is 0.0 wrt UTM grid)
41        .is_absolute    ??
42
43    """
44
45    ##
46    # @brief Instantiate an instance of class Geo_reference.
47    # @param zone The UTM zone.
48    # @param xllcorner X coord of origin of georef.
49    # @param yllcorner Y coord of origin of georef.
50    # @param datum ??
51    # @param projection The projection used (default UTM).
52    # @param units Units used in measuring distance (default m).
53    # @param false_easting ??
54    # @param false_northing ??
55    # @param NetCDFObject NetCDF file *handle* to write to.
56    # @param ASCIIFile ASCII text file *handle* to write to.
57    # @param read_title Title of the georeference text.
58    def __init__(self,
59                 zone=DEFAULT_ZONE,
60                 xllcorner=0.0,
61                 yllcorner=0.0,
62                 datum=DEFAULT_DATUM,
63                 projection=DEFAULT_PROJECTION,
64                 units=DEFAULT_UNITS,
65                 false_easting=DEFAULT_FALSE_EASTING,
66                 false_northing=DEFAULT_FALSE_NORTHING,
67                 NetCDFObject=None,
68                 ASCIIFile=None,
69                 read_title=None):
70        """
71        input:
72        NetCDFObject - a handle to the netCDF file to be written to
73        ASCIIFile - a handle to the text file
74        read_title - the title of the georeference text, if it was read in.
75         If the function that calls this has already read the title line,
76         it can't unread it, so this info has to be passed.
77         If you know of a way to unread this info, then tell us.
78
79         Note, the text file only saves a sub set of the info the
80         points file does.  Currently the info not written in text
81         must be the default info, since ANUGA assumes it isn't
82         changing.
83        """
84
85        if zone is None:
86            zone = DEFAULT_ZONE
87        self.false_easting = false_easting
88        self.false_northing = false_northing
89        self.datum = datum
90        self.projection = projection
91        self.zone = zone
92        self.units = units
93        self.xllcorner = xllcorner
94        self.yllcorner = yllcorner
95        #self.is_absolute = num.allclose([self.xllcorner, self.yllcorner], 0)
96
97        if NetCDFObject is not None:
98            self.read_NetCDF(NetCDFObject)
99
100        if ASCIIFile is not None:
101            self.read_ASCII(ASCIIFile, read_title=read_title)
102
103#    # Might be better to have this method instead of the following 3.
104#    def get_origin(self):
105#        return (self.zone, self.xllcorner, self.yllcorner)
106
107    ##
108    # @brief Get the X coordinate of the origin of this georef.
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################################################################################
239
240    def change_points_geo_ref(self, points, points_geo_ref=None):
241        """
242        Change the geo reference of a list or numeric array of points to
243        be this reference.(The reference used for this object)
244        If the points do not have a geo ref, assume 'absolute' values
245        """
246
247        is_list = False
248        if type(points) == types.ListType:
249            is_list = True
250
251        points = ensure_numeric(points, num.float)
252
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               '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               '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            if not points_geo_ref is None:
272                # Convert points to absolute coordinates
273                points[:,0] += points_geo_ref.xllcorner
274                points[:,1] += points_geo_ref.yllcorner
275
276            # Make points relative to primary geo reference
277            points[:,0] -= self.xllcorner
278            points[:,1] -= self.yllcorner
279
280        # if we got a list, return a list
281        if is_list:
282            points = points.tolist()
283
284        return points
285
286    ##
287    # @brief Is this georef absolute?
288    # @return True if ???
289    def is_absolute(self):
290        """Return True if xllcorner==yllcorner==0
291        indicating that points in question are absolute.
292        """
293
294        return num.allclose([self.xllcorner, self.yllcorner], 0)
295
296    ##
297    # @brief Convert points to absolute points in this georef instance.
298    # @param points
299    # @return
300    # @note
301    def get_absolute(self, points):
302        """Given a set of points geo referenced to this instance,
303        return the points as absolute values.
304        """
305
306        #if self.is_absolute:
307        #    return points
308
309        # remember if we got a list
310        is_list = False
311        if type(points) == types.ListType:
312            is_list = True
313
314        points = ensure_numeric(points, num.float)
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                #points = reshape(points, (1,2))
321
322        if not points.shape[1] == 2:
323            msg = ('Input must be an N x 2 array or list of (x,y) values. '
324                   'I got an %d x %d array' % points.shape)
325            raise ShapeError, msg
326
327        # Add geo ref to points
328        #if not self.is_absolute:
329        if not self.is_absolute():
330            points[:,0] += self.xllcorner
331            points[:,1] += self.yllcorner
332            #self.is_absolute = True
333
334        # if we got a list, return a list
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        is_list = False
353        if type(points) == types.ListType:
354            is_list = True
355
356        points = ensure_numeric(points, num.float)
357        if len(points.shape) == 1:
358            #One point has been passed
359            if not len(points) == 2:
360                msg = 'Single point must have two elements'
361                raise ShapeError, msg
362                #points = reshape(points, (1,2))
363
364        if not points.shape[1] == 2:
365            msg = ('Input must be an N x 2 array or list of (x,y) values. '
366                   'I got an %d x %d array' % points.shape)
367            raise ShapeError, msg
368
369        # Subtract geo ref from points
370        if not self.is_absolute():
371            points[:,0] -= self.xllcorner
372            points[:,1] -= self.yllcorner
373
374        # if we got a list, return a list
375        if is_list:
376            points = points.tolist()
377
378        return points
379
380    ##
381    # @brief ??
382    # @param other ??
383    def reconcile_zones(self, other):
384        if other is None:
385            other = Geo_reference()
386        if (self.zone == other.zone or
387            self.zone == DEFAULT_ZONE and
388            other.zone == DEFAULT_ZONE):
389            pass
390        elif self.zone == DEFAULT_ZONE:
391            self.zone = other.zone
392        elif other.zone == DEFAULT_ZONE:
393            other.zone = self.zone
394        else:
395            msg = ('Geospatial data must be in the same '
396                   'ZONE to allow reconciliation. I got zone %d and %d'
397                   % (self.zone, other.zone))
398            raise ANUGAError, msg
399
400    #def easting_northing2geo_reffed_point(self, x, y):
401    #    return [x-self.xllcorner, y - self.xllcorner]
402
403    #def easting_northing2geo_reffed_points(self, x, y):
404    #    return [x-self.xllcorner, y - self.xllcorner]
405
406    ##
407    # @brief Get origin of this geo_reference.
408    # @return (zone, xllcorner, yllcorner).
409    def get_origin(self):
410        return (self.zone, self.xllcorner, self.yllcorner)
411
412    ##
413    # @brief Get a string representation of this geo_reference instance.
414    def __repr__(self):
415        return ('(zone=%i easting=%f, northing=%f)'
416                % (self.zone, self.xllcorner, self.yllcorner))
417
418    ##
419    # @brief Compare two geo_reference instances.
420    # @param self This geo_reference instance.
421    # @param other Another geo_reference instance to compare against.
422    # @return 0 if instances have the same attributes, else 1.
423    # @note Attributes are: zone, xllcorner, yllcorner.
424    def __cmp__(self, other):
425        # FIXME (DSG) add a tolerence
426        if other is None:
427            return 1
428        cmp = 0
429        if not (self.xllcorner == self.xllcorner):
430            cmp = 1
431        if not (self.yllcorner == self.yllcorner):
432            cmp = 1
433        if not (self.zone == self.zone):
434            cmp = 1
435        return cmp
436
437
438##
439# @brief Write a geo_reference to a NetCDF file (usually SWW).
440# @param origin A georef instance or parameters to create a georef instance.
441# @param outfile Path to file to write.
442# @return A normalized geo_reference.
443def write_NetCDF_georeference(origin, outfile):
444    """Write georeference info to a netcdf file, usually sww.
445
446    The origin can be a georef instance or parameters for a geo_ref instance
447
448    outfile is the name of the file to be written to.
449    """
450
451    geo_ref = ensure_geo_reference(origin)
452    geo_ref.write_NetCDF(outfile)
453    return geo_ref
454
455
456##
457# @brief Convert an object to a georeference instance.
458# @param origin A georef instance or (zone, xllcorner, yllcorner)
459# @return A georef object, or None if 'origin' was None.
460def ensure_geo_reference(origin):
461    """
462    Given a list/tuple of zone, xllcorner and yllcorner of a geo-ref object,
463    return a geo ref object.
464
465    If the origin is None, return None, so calling this function doesn't
466    effect code logic
467    """
468
469    if isinstance(origin, Geo_reference):
470        geo_ref = origin
471    elif origin is None:
472        geo_ref = None
473    else:
474        geo_ref = apply(Geo_reference, origin)
475
476    return geo_ref
477
478
479#-----------------------------------------------------------------------
480
481if __name__ == "__main__":
482    pass
Note: See TracBrowser for help on using the repository browser.