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

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

Implemented run-time version of get_energy_through_cross_section and tested it.

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 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 = 'Points array must be two dimensional.\n'
207        msg += 'I got %d dimensions' %len(points.shape)
208        assert len(points.shape) == 2, msg
209
210        msg = 'Input must be an N x 2 array or list of (x,y) values. '
211        msg += 'I got an %d x %d array' %points.shape   
212        assert points.shape[1] == 2, msg               
213
214       
215        # FIXME (Ole): Could also check if zone, xllcorner, yllcorner
216        # are identical in the two geo refs.   
217        if points_geo_ref is not self:
218            # If georeferences are different
219       
220            if not points_geo_ref is None:
221                # Convert points to absolute coordinates
222                points[:,0] += points_geo_ref.xllcorner
223                points[:,1] += points_geo_ref.yllcorner
224       
225            # Make points relative to primary geo reference
226            points[:,0] -= self.xllcorner
227            points[:,1] -= self.yllcorner
228
229           
230        if is_list:
231            points = points.tolist()
232           
233        return points
234
235   
236    def is_absolute(self):
237        """Return True if xllcorner==yllcorner==0 indicating that points
238        in question are absolute.
239        """
240
241        return allclose([self.xllcorner, self.yllcorner], 0) 
242
243       
244   
245    def get_absolute(self, points):
246        """
247        Given a set of points geo referenced to this instance,
248        return the points as absolute values.
249        """
250
251        #if self.is_absolute():
252        #    return points
253        is_list = False
254        if type(points) == types.ListType:
255            is_list = True
256
257        points = ensure_numeric(points, Float)
258        if len(points.shape) == 1:
259            #One point has been passed
260            msg = 'Single point must have two elements'
261            if not len(points) == 2:
262                raise ShapeError, msg   
263                #points = reshape(points, (1,2))
264
265
266        msg = 'Input must be an N x 2 array or list of (x,y) values. '
267        msg += 'I got an %d x %d array' %points.shape   
268        if not points.shape[1] == 2:
269            raise ShapeError, msg   
270           
271       
272        # Add geo ref to points
273        if not self.is_absolute():
274            points[:,0] += self.xllcorner
275            points[:,1] += self.yllcorner
276
277       
278        if is_list:
279            points = points.tolist()
280             
281        return points
282
283
284    def get_relative(self, points):
285        """
286        Given a set of points in absolute UTM coordinates,
287        make them relative to this geo_reference instance,
288        return the points as relative values.
289
290        This is the inverse of get_absolute.
291        """
292
293        is_list = False
294        if type(points) == types.ListType:
295            is_list = True
296
297        points = ensure_numeric(points, Float)
298        if len(points.shape) == 1:
299            #One point has been passed
300            msg = 'Single point must have two elements'
301            if not len(points) == 2:
302                raise ShapeError, msg   
303                #points = reshape(points, (1,2))
304
305
306        msg = 'Input must be an N x 2 array or list of (x,y) values. '
307        msg += 'I got an %d x %d array' %points.shape   
308        if not points.shape[1] == 2:
309            raise ShapeError, msg   
310           
311       
312        # Subtract geo ref from points
313        if not self.is_absolute():
314            points[:,0] -= self.xllcorner
315            points[:,1] -= self.yllcorner
316
317       
318        if is_list:
319            points = points.tolist()
320             
321        return points
322
323
324
325    def reconcile_zones(self, other):
326        if other is None:
327            other = Geo_reference()
328        if self.zone == other.zone or\
329               self.zone == DEFAULT_ZONE and other.zone == DEFAULT_ZONE:
330            pass
331        elif self.zone == DEFAULT_ZONE:
332            self.zone = other.zone
333        elif other.zone == DEFAULT_ZONE:
334            other.zone = self.zone           
335        else:   
336            msg = 'Geospatial data must be in the same '+\
337                  'ZONE to allow reconciliation. I got zone %d and %d'\
338                  %(self.zone, other.zone)
339            raise ANUGAError, msg
340   
341    #def easting_northing2geo_reffed_point(self, x, y):
342    #    return [x-self.xllcorner, y - self.xllcorner]
343
344    #def easting_northing2geo_reffed_points(self, x, y):
345    #    return [x-self.xllcorner, y - self.xllcorner]
346
347    def get_origin(self):
348        return (self.zone, self.xllcorner, self.yllcorner)
349   
350    def __repr__(self):
351        return "(zone=%i easting=%f, northing=%f)" %(self.zone, self.xllcorner, self.yllcorner)
352   
353    def __cmp__(self,other):
354        # FIXME (DSG) add a tolerence
355        if other is None:
356            return 1
357        cmp = 0
358        if not (self.xllcorner == self.xllcorner):
359            cmp = 1
360        if not (self.yllcorner == self.yllcorner):
361            cmp = 1
362        if not (self.zone == self.zone):
363            cmp = 1
364        return cmp
365
366def write_NetCDF_georeference(origin, outfile):
367    """
368    Write georeferrence info to a netcdf file, usually sww.
369
370    The origin can be a georef instance or parrameters for a geo_ref instance
371
372    outfile is the name of the file to be written to.
373    """
374    geo_ref = ensure_geo_reference(origin)
375    geo_ref.write_NetCDF(outfile)
376    return geo_ref
377
378def ensure_geo_reference(origin):
379    """
380    Given a list/tuple of zone, xllcorner and yllcorner of a geo-ref object,
381    return a geo ref object.
382
383    If the origin is None, return None, so calling this function doesn't
384    effect code logic
385    """
386    if isinstance(origin, Geo_reference): 
387        geo_ref = origin
388    elif origin is None:
389        geo_ref = None
390    else:
391        geo_ref = apply(Geo_reference, origin)       
392    return geo_ref
393
394   
395#-----------------------------------------------------------------------
396
397if __name__ == "__main__":
398    pass
Note: See TracBrowser for help on using the repository browser.