source: anuga_core/source/anuga/coordinate_transforms/geo_reference.py @ 7062

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

Back-merge of changes to fix validation.

File size: 13.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
9
10import types, sys
11import copy
12from anuga.utilities.numerical_tools import ensure_numeric
13from anuga.utilities.anuga_exceptions import ANUGAError, TitleError, ParsingError, \
14     ShapeError
15
16import Numeric as num
17
18
19DEFAULT_ZONE = -1
20TITLE = '#geo reference' + "\n" #this title is referred to in the test format
21
22DEFAULT_PROJECTION = 'UTM'
23DEFAULT_DATUM = 'wgs84'
24DEFAULT_UNITS = 'm'
25DEFAULT_FALSE_EASTING = 500000
26DEFAULT_FALSE_NORTHING = 10000000 #Default for southern hemisphere
27
28class Geo_reference:
29    """
30    """
31
32    def __init__(self,
33                 zone = DEFAULT_ZONE,
34                 xllcorner = 0.0,
35                 yllcorner = 0.0,
36                 datum = DEFAULT_DATUM,
37                 projection = DEFAULT_PROJECTION,
38                 units = DEFAULT_UNITS,
39                 false_easting = DEFAULT_FALSE_EASTING,
40                 false_northing = DEFAULT_FALSE_NORTHING, 
41                 NetCDFObject=None,
42                 ASCIIFile=None,
43                 read_title=None):
44        """
45        input:
46        NetCDFObject - a handle to the netCDF file to be written to
47        ASCIIFile - a handle to the text file
48        read_title - the title of the georeference text, if it was read in.
49         If the function that calls this has already read the title line,
50         it can't unread it, so this info has to be passed.
51         If you know of a way to unread this info, then tell us.
52
53         Note, the text file only saves a sub set of the info the
54         points file does.  Currently the info not written in text
55         must be the default info, since ANUGA assumes it isn't
56         changing.
57         """
58        if zone is None:
59            zone = DEFAULT_ZONE
60        self.false_easting = int(false_easting)
61        self.false_northing = int(false_northing)
62        self.datum = datum
63        self.projection = projection
64        self.zone = int(zone)
65        self.units = units
66        self.xllcorner = float(xllcorner)
67        self.yllcorner = float(yllcorner)
68           
69        if NetCDFObject is not None:
70            self.read_NetCDF(NetCDFObject)
71 
72        if ASCIIFile is not None:
73            self.read_ASCII(ASCIIFile, read_title=read_title)
74           
75        # Set flag for absolute points (used by get_absolute)   
76        self.absolute = num.allclose([self.xllcorner, self.yllcorner], 0)
77           
78
79    def get_xllcorner(self):
80        return self.xllcorner
81   
82    def get_yllcorner(self):
83        return self.yllcorner
84   
85    def get_zone(self):
86        return self.zone
87       
88    def write_NetCDF(self, outfile):
89        outfile.xllcorner = self.xllcorner
90        outfile.yllcorner = self.yllcorner
91        outfile.zone = self.zone
92
93        outfile.false_easting = self.false_easting
94        outfile.false_northing = self.false_northing
95
96        outfile.datum = self.datum       
97        outfile.projection = self.projection
98        outfile.units = self.units
99
100    def read_NetCDF(self, infile):
101        self.xllcorner = float(infile.xllcorner[0])
102        self.yllcorner = float(infile.yllcorner[0])
103        self.zone = int(infile.zone[0])
104
105        try:
106            self.false_easting = int(infile.false_easting[0])
107            self.false_northing = int(infile.false_northing[0])
108
109            self.datum = infile.datum
110            self.projection = infile.projection
111            self.units = infile.units
112        except:
113            pass
114
115        if (self.false_easting != DEFAULT_FALSE_EASTING):
116            print "WARNING: False easting of %f specified." %self.false_easting
117            print "Default false easting is %f." %DEFAULT_FALSE_EASTING
118            print "ANUGA does not correct for differences in False Eastings."
119           
120        if (self.false_northing != DEFAULT_FALSE_NORTHING):
121            print "WARNING: False northing of %f specified." \
122                  %self.false_northing
123            print "Default false northing is %f." %DEFAULT_FALSE_NORTHING
124            print "ANUGA does not correct for differences in False Northings."
125           
126        if (self.datum.upper() != DEFAULT_DATUM.upper()):
127            print "WARNING: Datum of %s specified." \
128                  %self.datum
129            print "Default Datum is %s." %DEFAULT_DATUM
130            print "ANUGA does not correct for differences in datums."
131           
132        if (self.projection.upper() != DEFAULT_PROJECTION.upper()):
133            print "WARNING: Projection of %s specified." \
134                  %self.projection
135            print "Default Projection is %s." %DEFAULT_PROJECTION
136            print "ANUGA does not correct for differences in Projection."
137           
138        if (self.units.upper() != DEFAULT_UNITS.upper()):
139            print "WARNING: Units of %s specified." \
140                  %self.units
141            print "Default units is %s." %DEFAULT_UNITS
142            print "ANUGA does not correct for differences in units."
143       
144    ### ASCII files with geo-refs are currently not used   
145    def write_ASCII(self, fd):
146        fd.write(TITLE)
147        fd.write(str(self.zone) + "\n") 
148        fd.write(str(self.xllcorner) + "\n") 
149        fd.write(str(self.yllcorner) + "\n")
150
151    def read_ASCII(self, fd,read_title=None):
152        try:
153            if read_title == None:
154                read_title = fd.readline() # remove the title line
155            if read_title[0:2].upper() != TITLE[0:2].upper():
156                msg = 'File error.  Expecting line: %s.  Got this line: %s' \
157                      %(TITLE, read_title)
158                raise TitleError, msg
159            self.zone = int(fd.readline())
160            self.xllcorner = float(fd.readline())
161            self.yllcorner = float(fd.readline())
162        except SyntaxError:
163                msg = 'File error.  Got syntax error while parsing geo reference' 
164                raise ParsingError, msg
165           
166        # Fix some assertion failures
167        if(type(self.zone) == num.ArrayType and self.zone.shape == ()):
168            self.zone = self.zone[0]
169        if type(self.xllcorner) == num.ArrayType and self.xllcorner.shape == ():
170            self.xllcorner = self.xllcorner[0]
171        if type(self.yllcorner) == num.ArrayType and self.yllcorner.shape == ():
172            self.yllcorner = self.yllcorner[0]
173   
174        assert (type(self.xllcorner) == types.FloatType)
175        assert (type(self.yllcorner) == types.FloatType)
176        assert (type(self.zone) == types.IntType)
177       
178       
179    def change_points_geo_ref(self, points,
180                              points_geo_ref=None):
181        """
182        Change the geo reference of a list or Numeric array of points to
183        be this reference.(The reference used for this object)
184        If the points do not have a geo ref, assume 'absolute' values
185        """
186        import copy
187       
188        is_list = False
189        if type(points) == types.ListType:
190            is_list = True
191
192        points = ensure_numeric(points, num.Float)
193       
194        if len(points.shape) == 1:
195            #One point has been passed
196            msg = 'Single point must have two elements'
197            assert len(points) == 2, msg
198            points = num.reshape(points, (1,2))
199
200        msg = 'Points array must be two dimensional.\n'
201        msg += 'I got %d dimensions' %len(points.shape)
202        assert len(points.shape) == 2, msg
203
204        msg = 'Input must be an N x 2 array or list of (x,y) values. '
205        msg += 'I got an %d x %d array' %points.shape   
206        assert points.shape[1] == 2, msg               
207
208       
209        # FIXME (Ole): Could also check if zone, xllcorner, yllcorner
210        # are identical in the two geo refs.   
211        if points_geo_ref is not self:
212            # If georeferences are different
213       
214            points = copy.copy(points) # Don't destroy input                   
215            if not points_geo_ref is None:
216                # Convert points to absolute coordinates
217                points[:,0] += points_geo_ref.xllcorner
218                points[:,1] += points_geo_ref.yllcorner
219       
220            # Make points relative to primary geo reference
221            points[:,0] -= self.xllcorner
222            points[:,1] -= self.yllcorner
223
224           
225        if is_list:
226            points = points.tolist()
227           
228        return points
229
230   
231    def is_absolute(self):
232        """Return True if xllcorner==yllcorner==0 indicating that points
233        in question are absolute.
234        """
235       
236        # FIXME(Ole): It is unfortunate that decision about whether points
237        # are absolute or not lies with the georeference object. Ross pointed this out.
238        # Moreover, this little function is responsible for a large fraction of the time
239        # using in data fitting (something in like 40 - 50%.
240        # This was due to the repeated calls to allclose.
241        # With the flag method fitting is much faster (18 Mar 2009).
242
243       
244        # FIXME(Ole): HACK to be able to reuse data already cached (18 Mar 2009).
245        # Remove at some point
246        if not hasattr(self, 'absolute'):
247            self.absolute = num.allclose([self.xllcorner, self.yllcorner], 0)
248
249           
250        # Return absolute flag   
251        return self.absolute
252       
253   
254    def get_absolute(self, points):
255        """
256        Given a set of points geo referenced to this instance,
257        return the points as absolute values.
258        """
259
260        is_list = False
261        if type(points) == types.ListType:
262            is_list = True
263
264        points = ensure_numeric(points, num.Float)
265        if len(points.shape) == 1:
266            # One point has been passed
267            msg = 'Single point must have two elements'
268            if not len(points) == 2:
269                raise ShapeError, msg   
270
271
272        msg = 'Input must be an N x 2 array or list of (x,y) values. '
273        msg += 'I got an %d x %d array' %points.shape   
274        if not points.shape[1] == 2:
275            raise ShapeError, msg   
276           
277       
278        # Add geo ref to points
279        if not self.is_absolute():
280            points = copy.copy(points) # Don't destroy input                   
281            points[:,0] += self.xllcorner
282            points[:,1] += self.yllcorner
283
284       
285        if is_list:
286            points = points.tolist()
287             
288        return points
289
290
291    def get_relative(self, points):
292        """
293        Given a set of points in absolute UTM coordinates,
294        make them relative to this geo_reference instance,
295        return the points as relative values.
296
297        This is the inverse of get_absolute.
298        """
299
300        is_list = False
301        if type(points) == types.ListType:
302            is_list = True
303
304        points = ensure_numeric(points, num.Float)
305        if len(points.shape) == 1:
306            #One point has been passed
307            msg = 'Single point must have two elements'
308            if not len(points) == 2:
309                raise ShapeError, msg   
310                #points = reshape(points, (1,2))
311
312
313        msg = 'Input must be an N x 2 array or list of (x,y) values. '
314        msg += 'I got an %d x %d array' %points.shape   
315        if not points.shape[1] == 2:
316            raise ShapeError, msg   
317           
318       
319        # Subtract geo ref from points
320        if not self.is_absolute():
321            points = copy.copy(points) # Don't destroy input                           
322            points[:,0] -= self.xllcorner
323            points[:,1] -= self.yllcorner
324
325       
326        if is_list:
327            points = points.tolist()
328             
329        return points
330
331
332
333    def reconcile_zones(self, other):
334        if other is None:
335            other = Geo_reference()
336        if self.zone == other.zone or\
337               self.zone == DEFAULT_ZONE and other.zone == DEFAULT_ZONE:
338            pass
339        elif self.zone == DEFAULT_ZONE:
340            self.zone = other.zone
341        elif other.zone == DEFAULT_ZONE:
342            other.zone = self.zone           
343        else:   
344            msg = 'Geospatial data must be in the same '+\
345                  'ZONE to allow reconciliation. I got zone %d and %d'\
346                  %(self.zone, other.zone)
347            raise ANUGAError, msg
348   
349    #def easting_northing2geo_reffed_point(self, x, y):
350    #    return [x-self.xllcorner, y - self.xllcorner]
351
352    #def easting_northing2geo_reffed_points(self, x, y):
353    #    return [x-self.xllcorner, y - self.xllcorner]
354
355    def get_origin(self):
356        return (self.zone, self.xllcorner, self.yllcorner)
357   
358    def __repr__(self):
359        return "(zone=%i easting=%f, northing=%f)" %(self.zone, self.xllcorner, self.yllcorner)
360   
361    def __cmp__(self,other):
362        # FIXME (DSG) add a tolerence
363        if other is None:
364            return 1
365        cmp = 0
366        if not (self.xllcorner == self.xllcorner):
367            cmp = 1
368        if not (self.yllcorner == self.yllcorner):
369            cmp = 1
370        if not (self.zone == self.zone):
371            cmp = 1
372        return cmp
373
374def write_NetCDF_georeference(origin, outfile):
375    """
376    Write georeferrence info to a netcdf file, usually sww.
377
378    The origin can be a georef instance or parrameters for a geo_ref instance
379
380    outfile is the name of the file to be written to.
381    """
382    geo_ref = ensure_geo_reference(origin)
383    geo_ref.write_NetCDF(outfile)
384    return geo_ref
385
386def ensure_geo_reference(origin):
387    """
388    Given a list/tuple of zone, xllcorner and yllcorner of a geo-ref object,
389    return a geo ref object.
390
391    If the origin is None, return None, so calling this function doesn't
392    effect code logic
393    """
394    if isinstance(origin, Geo_reference): 
395        geo_ref = origin
396    elif origin is None:
397        geo_ref = None
398    else:
399        geo_ref = apply(Geo_reference, origin)       
400    return geo_ref
401
402   
403#-----------------------------------------------------------------------
404
405if __name__ == "__main__":
406    pass
Note: See TracBrowser for help on using the repository browser.