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

Last change on this file since 4387 was 4387, checked in by duncan, 17 years ago

small modifications to geo-reffing
File size: 10.1 KB
Line 
1"""
2"""
3
4
5#FIXME: Ensure that all attributes of a georef are treated everywhere
6#and unit test
7
8import types, sys
9from Numeric import array, Float, ArrayType, reshape, allclose
10from anuga.utilities.numerical_tools import ensure_numeric
11from anuga.utilities.anuga_exceptions import ANUGAError, TitleError, ParsingError, \
12     ShapeError
13
14
15DEFAULT_ZONE = -1
16TITLE = '#geo reference' + "\n" #this title is referred to in the .xya format
17
18DEFAULT_PROJECTION = 'UTM'
19DEFAULT_DATUM = 'wgs84'
20DEFAULT_UNITS = 'm'
21DEFAULT_FALSE_EASTING = 500000
22DEFAULT_FALSE_NORTHING = 10000000 #Default for southern hemisphere
23
24class Geo_reference:
25    """
26    """
27
28    def __init__(self,
29                 zone = DEFAULT_ZONE,
30                 xllcorner = 0.0,
31                 yllcorner = 0.0,
32                 datum = DEFAULT_DATUM,
33                 projection = DEFAULT_PROJECTION,
34                 units = DEFAULT_UNITS,
35                 false_easting = DEFAULT_FALSE_EASTING,
36                 false_northing = DEFAULT_FALSE_NORTHING, 
37                 NetCDFObject=None,
38                 ASCIIFile=None,
39                 read_title=None):
40        """
41        input:
42        NetCDFObject - a handle to the netCDF file to be written to
43        ASCIIFile - a handle to the text file
44        read_title - the title of the georeference text, if it was read in.
45         If the function that calls this has already read the title line,
46         it can't unread it, so this info has to be passed.
47         If you know of a way to unread this info, then tell us.
48
49         Note, the text file only saves a sub set of the info the
50         points file does.  Currently the info not written in text
51         must be the default info, since ANUGA assumes it isn't
52         changing.
53         """
54
55        self.false_easting = false_easting
56        self.false_northing = false_northing       
57        self.datum = datum
58        self.projection = projection
59        self.zone = zone       
60        self.units = units
61        self.xllcorner = xllcorner
62        self.yllcorner = yllcorner       
63           
64        if NetCDFObject is not None:
65            self.read_NetCDF(NetCDFObject)
66 
67        if ASCIIFile is not None:
68            self.read_ASCII(ASCIIFile, read_title=read_title)
69
70    def get_xllcorner(self):
71        return self.xllcorner
72   
73    def get_yllcorner(self):
74        return self.yllcorner
75   
76    def get_zone(self):
77        return self.zone
78       
79    def write_NetCDF(self, outfile):
80        outfile.xllcorner = self.xllcorner
81        outfile.yllcorner = self.yllcorner
82        outfile.zone = self.zone
83
84        outfile.false_easting = self.false_easting
85        outfile.false_northing = self.false_northing
86
87        outfile.datum = self.datum       
88        outfile.projection = self.projection
89        outfile.units = self.units
90
91    def read_NetCDF(self, infile):
92        self.xllcorner = infile.xllcorner[0]
93        self.yllcorner = infile.yllcorner[0] 
94        self.zone = infile.zone[0]
95
96        self.false_easting = infile.false_easting[0]
97        self.false_northing = infile.false_northing[0]
98       
99        self.datum = infile.datum       
100        self.projection = infile.projection
101        self.units = infile.units
102       
103        # Fix some assertion failures
104        if type(self.zone) == ArrayType and self.zone.shape == ():
105            self.zone = self.zone[0]
106        if type(self.xllcorner) == ArrayType and self.xllcorner.shape == ():
107            self.xllcorner = self.xllcorner[0]
108        if type(self.yllcorner) == ArrayType and self.yllcorner.shape == ():
109            self.yllcorner = self.yllcorner[0]
110
111        assert (type(self.xllcorner) == types.FloatType or\
112                type(self.xllcorner) == types.IntType)
113        assert (type(self.yllcorner) == types.FloatType or\
114                type(self.yllcorner) == types.IntType)
115        assert (type(self.zone) == types.IntType)
116       
117        assert (self.false_easting == DEFAULT_FALSE_EASTING)
118        assert (self.false_northing == DEFAULT_FALSE_NORTHING)
119
120        assert(self.datum == DEFAULT_DATUM)
121        assert(self.projection == DEFAULT_PROJECTION)
122        assert (self.units == DEFAULT_UNITS)
123       
124       
125    def write_ASCII(self, fd):
126        fd.write(TITLE)
127        fd.write(str(self.zone) + "\n") 
128        fd.write(str(self.xllcorner) + "\n") 
129        fd.write(str(self.yllcorner) + "\n")
130
131    def read_ASCII(self, fd,read_title=None):
132        try:
133            if read_title == None:
134                read_title = fd.readline() # remove the title line
135            if read_title[0:2].upper() != TITLE[0:2].upper():
136                msg = 'File error.  Expecting line: %s.  Got this line: %s' \
137                      %(TITLE, read_title)
138                raise TitleError, msg
139            self.zone = int(fd.readline())
140            self.xllcorner = float(fd.readline())
141            self.yllcorner = float(fd.readline())
142        except SyntaxError:
143                msg = 'File error.  Got syntax error while parsing geo reference' 
144                raise ParsingError, msg
145           
146        # Fix some assertion failures
147        if(type(self.zone) == ArrayType and self.zone.shape == ()):
148            self.zone = self.zone[0]
149        if type(self.xllcorner) == ArrayType and self.xllcorner.shape == ():
150            self.xllcorner = self.xllcorner[0]
151        if type(self.yllcorner) == ArrayType and self.yllcorner.shape == ():
152            self.yllcorner = self.yllcorner[0]
153   
154        assert (type(self.xllcorner) == types.FloatType)
155        assert (type(self.yllcorner) == types.FloatType)
156        assert (type(self.zone) == types.IntType)
157       
158       
159    def change_points_geo_ref(self, points,
160                              points_geo_ref=None):
161        """
162        Change the geo reference of a list or Numeric array of points to
163        be this reference.(The reference used for this object)
164        If the points do not have a geo ref, assume 'absolute' values
165        """
166
167        is_list = False
168        if type(points) == types.ListType:
169            is_list = True
170
171        points = ensure_numeric(points, Float)
172       
173        if len(points.shape) == 1:
174            #One point has been passed
175            msg = 'Single point must have two elements'
176            assert len(points) == 2, msg
177            points = reshape(points, (1,2))
178
179        msg = 'Input must be an N x 2 array or list of (x,y) values. '
180        msg += 'I got an %d x %d array' %points.shape   
181        assert points.shape[1] == 2, msg               
182
183           
184        if points_geo_ref is not self:
185            #add point geo ref to points
186            if not points_geo_ref is None:
187                points[:,0] += points_geo_ref.xllcorner
188                points[:,1] += points_geo_ref.yllcorner
189       
190            #subtract primary geo ref from points
191            points[:,0] -= self.xllcorner
192            points[:,1] -= self.yllcorner
193
194           
195        if is_list:
196            points = points.tolist()
197           
198        return points
199
200   
201    def is_absolute(self):
202        """Return True if xllcorner==yllcorner==0 indicating that points
203        in question are absolute.
204        """
205
206        return allclose([self.xllcorner, self.yllcorner], 0) 
207
208       
209   
210    def get_absolute(self, points):
211        """
212        Given a set of points geo referenced to this instance,
213        return the points as absolute values.
214        """
215
216        #if self.is_absolute():
217        #    return points
218       
219
220        is_list = False
221        if type(points) == types.ListType:
222            is_list = True
223
224        points = ensure_numeric(points, Float)
225        if len(points.shape) == 1:
226            #One point has been passed
227            msg = 'Single point must have two elements'
228            if not len(points) == 2:
229                raise ShapeError, msg   
230                #points = reshape(points, (1,2))
231
232
233        msg = 'Input must be an N x 2 array or list of (x,y) values. '
234        msg += 'I got an %d x %d array' %points.shape   
235        if not points.shape[1] == 2:
236            raise ShapeError, msg   
237           
238       
239        # Add geo ref to points
240        if not self.is_absolute():
241            points[:,0] += self.xllcorner
242            points[:,1] += self.yllcorner
243
244       
245        if is_list:
246            points = points.tolist()
247             
248        return points
249
250
251    def reconcile_zones(self, other):
252        if other is None:
253            other = Geo_reference()
254        if self.zone == other.zone or\
255               self.zone == DEFAULT_ZONE and other.zone == DEFAULT_ZONE:
256            pass
257        elif self.zone == DEFAULT_ZONE:
258            self.zone = other.zone
259        elif other.zone == DEFAULT_ZONE:
260            other.zone = self.zone           
261        else:   
262            msg = 'Both geospatial_data objects must be in the same '+\
263                  'ZONE to allow reconciliation. I got zone %d and %d'\
264                  %(self.zone, other.zone)
265            raise ANUGAError, msg
266   
267    #def easting_northing2geo_reffed_point(self, x, y):
268    #    return [x-self.xllcorner, y - self.xllcorner]
269
270    #def easting_northing2geo_reffed_points(self, x, y):
271    #    return [x-self.xllcorner, y - self.xllcorner]
272
273    def get_origin(self):
274        return (self.zone, self.xllcorner, self.yllcorner)
275   
276    def __repr__(self):
277        return "(zone=%i easting=%f, northing=%f)" %(self.zone, self.xllcorner, self.yllcorner)
278   
279    def __cmp__(self,other):
280        # FIXME (DSG) add a tolerence
281        if other is None:
282            return 1
283        cmp = 0
284        if not (self.xllcorner == self.xllcorner):
285            cmp = 1
286        if not (self.yllcorner == self.yllcorner):
287            cmp = 1
288        if not (self.zone == self.zone):
289            cmp = 1
290        return cmp
291
292def write_NetCDF_georeference(origin, outfile):
293    """
294    Write georeferrence info to a netcdf file, usually sww.
295
296    The origin can be a georef instance or parrameters for a geo_ref instance
297
298    outfile is the name of the file to be written to.
299    """
300    if isinstance(origin, Geo_reference): 
301        geo_ref = origin
302    else:
303        geo_ref = apply(Geo_reference,origin)
304    geo_ref.write_NetCDF(outfile)
305    return geo_ref
306   
307#-----------------------------------------------------------------------
308
309if __name__ == "__main__":
310    pass
Note: See TracBrowser for help on using the repository browser.