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

Last change on this file since 4839 was 4839, checked in by duncan, 16 years ago

speeding up slow functions in general mesh. Fix for ticket 201. Also checking in GA validation test for profiling fitting
File size: 11.8 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 reconcile_zones(self, other):
277        if other is None:
278            other = Geo_reference()
279        if self.zone == other.zone or\
280               self.zone == DEFAULT_ZONE and other.zone == DEFAULT_ZONE:
281            pass
282        elif self.zone == DEFAULT_ZONE:
283            self.zone = other.zone
284        elif other.zone == DEFAULT_ZONE:
285            other.zone = self.zone           
286        else:   
287            msg = 'Geospatial data must be in the same '+\
288                  'ZONE to allow reconciliation. I got zone %d and %d'\
289                  %(self.zone, other.zone)
290            raise ANUGAError, msg
291   
292    #def easting_northing2geo_reffed_point(self, x, y):
293    #    return [x-self.xllcorner, y - self.xllcorner]
294
295    #def easting_northing2geo_reffed_points(self, x, y):
296    #    return [x-self.xllcorner, y - self.xllcorner]
297
298    def get_origin(self):
299        return (self.zone, self.xllcorner, self.yllcorner)
300   
301    def __repr__(self):
302        return "(zone=%i easting=%f, northing=%f)" %(self.zone, self.xllcorner, self.yllcorner)
303   
304    def __cmp__(self,other):
305        # FIXME (DSG) add a tolerence
306        if other is None:
307            return 1
308        cmp = 0
309        if not (self.xllcorner == self.xllcorner):
310            cmp = 1
311        if not (self.yllcorner == self.yllcorner):
312            cmp = 1
313        if not (self.zone == self.zone):
314            cmp = 1
315        return cmp
316
317def write_NetCDF_georeference(origin, outfile):
318    """
319    Write georeferrence info to a netcdf file, usually sww.
320
321    The origin can be a georef instance or parrameters for a geo_ref instance
322
323    outfile is the name of the file to be written to.
324    """
325    geo_ref = ensure_geo_reference(origin)
326    geo_ref.write_NetCDF(outfile)
327    return geo_ref
328
329def ensure_geo_reference(origin):
330    """
331    Given a list/tuple of zone, xllcorner and yllcorner of a geo-ref object,
332    return a geo ref object.
333
334    If the origin is None, return None, so calling this function doesn't
335    effect code logic
336    """
337    if isinstance(origin, Geo_reference): 
338        geo_ref = origin
339    elif origin is None:
340        geo_ref = None
341    else:
342        geo_ref = apply(Geo_reference,origin)       
343    return geo_ref
344
345   
346#-----------------------------------------------------------------------
347
348if __name__ == "__main__":
349    pass
Note: See TracBrowser for help on using the repository browser.