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

Last change on this file since 6149 was 6149, checked in by rwilson, 16 years ago

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