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

Last change on this file since 5421 was 5288, checked in by ole, 16 years ago

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