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

Last change on this file since 6438 was 6438, checked in by ole, 15 years ago

Fixed bug in ensure_numeric where input was always copied irrespective of type.

Then fixed snowball bugs by copying input to ensure_numeric where appropriate.

Thanks Ross

File size: 13.4 KB
RevLine 
[2488]1"""
[4451]2
[2488]3"""
4
5
6#FIXME: Ensure that all attributes of a georef are treated everywhere
7#and unit test
8
[6438]9
[2593]10import types, sys
[6438]11import copy
[3514]12from anuga.utilities.numerical_tools import ensure_numeric
13from anuga.utilities.anuga_exceptions import ANUGAError, TitleError, ParsingError, \
[3134]14     ShapeError
[2488]15
[6149]16import Numeric as num
[2941]17
[6149]18
[2602]19DEFAULT_ZONE = -1
[4663]20TITLE = '#geo reference' + "\n" #this title is referred to in the test format
[2488]21
[4387]22DEFAULT_PROJECTION = 'UTM'
23DEFAULT_DATUM = 'wgs84'
24DEFAULT_UNITS = 'm'
25DEFAULT_FALSE_EASTING = 500000
26DEFAULT_FALSE_NORTHING = 10000000 #Default for southern hemisphere
27
[2488]28class Geo_reference:
29    """
30    """
31
32    def __init__(self,
33                 zone = DEFAULT_ZONE,
34                 xllcorner = 0.0,
35                 yllcorner = 0.0,
[4387]36                 datum = DEFAULT_DATUM,
37                 projection = DEFAULT_PROJECTION,
38                 units = DEFAULT_UNITS,
39                 false_easting = DEFAULT_FALSE_EASTING,
40                 false_northing = DEFAULT_FALSE_NORTHING, 
[2488]41                 NetCDFObject=None,
42                 ASCIIFile=None,
43                 read_title=None):
44        """
45        input:
46        NetCDFObject - a handle to the netCDF file to be written to
47        ASCIIFile - a handle to the text file
48        read_title - the title of the georeference text, if it was read in.
[2941]49         If the function that calls this has already read the title line,
50         it can't unread it, so this info has to be passed.
51         If you know of a way to unread this info, then tell us.
[2488]52
[4387]53         Note, the text file only saves a sub set of the info the
54         points file does.  Currently the info not written in text
55         must be the default info, since ANUGA assumes it isn't
56         changing.
57         """
[4451]58        if zone is None:
59            zone = DEFAULT_ZONE
[2488]60        self.false_easting = false_easting
61        self.false_northing = false_northing       
62        self.datum = datum
63        self.projection = projection
64        self.zone = zone       
65        self.units = units
66        self.xllcorner = xllcorner
67        self.yllcorner = yllcorner       
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
[4387]101       
[2488]102        # Fix some assertion failures
[6149]103        if type(self.zone) == num.ArrayType and self.zone.shape == ():
[2488]104            self.zone = self.zone[0]
[6149]105        if type(self.xllcorner) == num.ArrayType and self.xllcorner.shape == ():
[2488]106            self.xllcorner = self.xllcorner[0]
[6149]107        if type(self.yllcorner) == num.ArrayType and self.yllcorner.shape == ():
[2488]108            self.yllcorner = self.yllcorner[0]
109
110        assert (type(self.xllcorner) == types.FloatType or\
111                type(self.xllcorner) == types.IntType)
112        assert (type(self.yllcorner) == types.FloatType or\
113                type(self.yllcorner) == types.IntType)
114        assert (type(self.zone) == types.IntType)
[4387]115       
[4389]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
[4451]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."
[2488]153       
[4663]154    ### ASCII files with geo-refs are currently not used   
[2488]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):
[2942]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
[2488]175           
176        # Fix some assertion failures
[6149]177        if(type(self.zone) == num.ArrayType and self.zone.shape == ()):
[2488]178            self.zone = self.zone[0]
[6149]179        if type(self.xllcorner) == num.ArrayType and self.xllcorner.shape == ():
[2488]180            self.xllcorner = self.xllcorner[0]
[6149]181        if type(self.yllcorner) == num.ArrayType and self.yllcorner.shape == ():
[2488]182            self.yllcorner = self.yllcorner[0]
183   
184        assert (type(self.xllcorner) == types.FloatType)
185        assert (type(self.yllcorner) == types.FloatType)
186        assert (type(self.zone) == types.IntType)
187       
188       
189    def change_points_geo_ref(self, points,
190                              points_geo_ref=None):
191        """
192        Change the geo reference of a list or Numeric array of points to
193        be this reference.(The reference used for this object)
194        If the points do not have a geo ref, assume 'absolute' values
195        """
[6438]196        import copy
197       
[2488]198        is_list = False
199        if type(points) == types.ListType:
200            is_list = True
[2594]201
[6438]202        points = ensure_numeric(copy.copy(points), num.Float)
[2595]203       
[2594]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
[6149]208            points = num.reshape(points, (1,2))
[2594]209
[5665]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
[3071]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
[5736]218       
219        # FIXME (Ole): Could also check if zone, xllcorner, yllcorner
220        # are identical in the two geo refs.   
[2594]221        if points_geo_ref is not self:
[5736]222            # If georeferences are different
223       
[2594]224            if not points_geo_ref is None:
[5736]225                # Convert points to absolute coordinates
[2594]226                points[:,0] += points_geo_ref.xllcorner
227                points[:,1] += points_geo_ref.yllcorner
[2488]228       
[5736]229            # Make points relative to primary geo reference
[2594]230            points[:,0] -= self.xllcorner
231            points[:,1] -= self.yllcorner
232
[2488]233           
234        if is_list:
235            points = points.tolist()
[2594]236           
[2488]237        return points
[3129]238
[2488]239   
[3129]240    def is_absolute(self):
241        """Return True if xllcorner==yllcorner==0 indicating that points
242        in question are absolute.
243        """
244
[6149]245        return num.allclose([self.xllcorner, self.yllcorner], 0) 
[3129]246
247       
[2594]248   
[2488]249    def get_absolute(self, points):
250        """
251        Given a set of points geo referenced to this instance,
252        return the points as absolute values.
253        """
[3071]254
[3134]255        #if self.is_absolute():
256        #    return points
[2488]257        is_list = False
258        if type(points) == types.ListType:
259            is_list = True
[3071]260
[6438]261        points = ensure_numeric(copy.copy(points), num.Float)
[3071]262        if len(points.shape) == 1:
263            #One point has been passed
264            msg = 'Single point must have two elements'
[3134]265            if not len(points) == 2:
266                raise ShapeError, msg   
267                #points = reshape(points, (1,2))
[3071]268
269
270        msg = 'Input must be an N x 2 array or list of (x,y) values. '
271        msg += 'I got an %d x %d array' %points.shape   
[3134]272        if not points.shape[1] == 2:
273            raise ShapeError, msg   
[2488]274           
[3071]275       
[3134]276        # Add geo ref to points
277        if not self.is_absolute():
278            points[:,0] += self.xllcorner
279            points[:,1] += self.yllcorner
[3071]280
281       
[2488]282        if is_list:
283            points = points.tolist()
284             
285        return points
286
[2594]287
[5288]288    def get_relative(self, points):
289        """
290        Given a set of points in absolute UTM coordinates,
291        make them relative to this geo_reference instance,
292        return the points as relative values.
293
294        This is the inverse of get_absolute.
295        """
296
297        is_list = False
298        if type(points) == types.ListType:
299            is_list = True
300
[6149]301        points = ensure_numeric(points, num.Float)
[5288]302        if len(points.shape) == 1:
303            #One point has been passed
304            msg = 'Single point must have two elements'
305            if not len(points) == 2:
306                raise ShapeError, msg   
307                #points = reshape(points, (1,2))
308
309
310        msg = 'Input must be an N x 2 array or list of (x,y) values. '
311        msg += 'I got an %d x %d array' %points.shape   
312        if not points.shape[1] == 2:
313            raise ShapeError, msg   
314           
315       
316        # Subtract geo ref from points
317        if not self.is_absolute():
318            points[:,0] -= self.xllcorner
319            points[:,1] -= self.yllcorner
320
321       
322        if is_list:
323            points = points.tolist()
324             
325        return points
326
327
328
[2594]329    def reconcile_zones(self, other):
[4180]330        if other is None:
331            other = Geo_reference()
[3071]332        if self.zone == other.zone or\
333               self.zone == DEFAULT_ZONE and other.zone == DEFAULT_ZONE:
[2594]334            pass
335        elif self.zone == DEFAULT_ZONE:
336            self.zone = other.zone
337        elif other.zone == DEFAULT_ZONE:
338            other.zone = self.zone           
339        else:   
[4470]340            msg = 'Geospatial data must be in the same '+\
[3614]341                  'ZONE to allow reconciliation. I got zone %d and %d'\
342                  %(self.zone, other.zone)
[2608]343            raise ANUGAError, msg
[2488]344   
345    #def easting_northing2geo_reffed_point(self, x, y):
346    #    return [x-self.xllcorner, y - self.xllcorner]
347
348    #def easting_northing2geo_reffed_points(self, x, y):
349    #    return [x-self.xllcorner, y - self.xllcorner]
350
351    def get_origin(self):
352        return (self.zone, self.xllcorner, self.yllcorner)
353   
354    def __repr__(self):
355        return "(zone=%i easting=%f, northing=%f)" %(self.zone, self.xllcorner, self.yllcorner)
356   
357    def __cmp__(self,other):
358        # FIXME (DSG) add a tolerence
359        if other is None:
360            return 1
361        cmp = 0
362        if not (self.xllcorner == self.xllcorner):
363            cmp = 1
364        if not (self.yllcorner == self.yllcorner):
365            cmp = 1
366        if not (self.zone == self.zone):
367            cmp = 1
368        return cmp
[4387]369
370def write_NetCDF_georeference(origin, outfile):
371    """
372    Write georeferrence info to a netcdf file, usually sww.
373
374    The origin can be a georef instance or parrameters for a geo_ref instance
375
376    outfile is the name of the file to be written to.
377    """
[4451]378    geo_ref = ensure_geo_reference(origin)
379    geo_ref.write_NetCDF(outfile)
380    return geo_ref
381
382def ensure_geo_reference(origin):
383    """
384    Given a list/tuple of zone, xllcorner and yllcorner of a geo-ref object,
385    return a geo ref object.
386
387    If the origin is None, return None, so calling this function doesn't
388    effect code logic
389    """
[4387]390    if isinstance(origin, Geo_reference): 
391        geo_ref = origin
[4451]392    elif origin is None:
393        geo_ref = None
[4387]394    else:
[5730]395        geo_ref = apply(Geo_reference, origin)       
[4387]396    return geo_ref
[4451]397
[4387]398   
[2488]399#-----------------------------------------------------------------------
400
401if __name__ == "__main__":
402    pass
Note: See TracBrowser for help on using the repository browser.