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

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

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