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

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

numpy changes

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