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

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

backwards compatibility

File size: 10.2 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       
97        # Fix some assertion failures
98        if type(self.zone) == ArrayType and self.zone.shape == ():
99            self.zone = self.zone[0]
100        if type(self.xllcorner) == ArrayType and self.xllcorner.shape == ():
101            self.xllcorner = self.xllcorner[0]
102        if type(self.yllcorner) == ArrayType and self.yllcorner.shape == ():
103            self.yllcorner = self.yllcorner[0]
104
105        assert (type(self.xllcorner) == types.FloatType or\
106                type(self.xllcorner) == types.IntType)
107        assert (type(self.yllcorner) == types.FloatType or\
108                type(self.yllcorner) == types.IntType)
109        assert (type(self.zone) == types.IntType)
110       
111        try:
112            self.false_easting = infile.false_easting[0]
113            self.false_northing = infile.false_northing[0]
114       
115            self.datum = infile.datum       
116            self.projection = infile.projection
117            self.units = infile.units
118        except:
119            pass
120        assert (self.false_easting == DEFAULT_FALSE_EASTING)
121        assert (self.false_northing == DEFAULT_FALSE_NORTHING)
122
123        assert(self.datum == DEFAULT_DATUM)
124        assert(self.projection == DEFAULT_PROJECTION)
125        assert (self.units == DEFAULT_UNITS)
126       
127       
128    def write_ASCII(self, fd):
129        fd.write(TITLE)
130        fd.write(str(self.zone) + "\n") 
131        fd.write(str(self.xllcorner) + "\n") 
132        fd.write(str(self.yllcorner) + "\n")
133
134    def read_ASCII(self, fd,read_title=None):
135        try:
136            if read_title == None:
137                read_title = fd.readline() # remove the title line
138            if read_title[0:2].upper() != TITLE[0:2].upper():
139                msg = 'File error.  Expecting line: %s.  Got this line: %s' \
140                      %(TITLE, read_title)
141                raise TitleError, msg
142            self.zone = int(fd.readline())
143            self.xllcorner = float(fd.readline())
144            self.yllcorner = float(fd.readline())
145        except SyntaxError:
146                msg = 'File error.  Got syntax error while parsing geo reference' 
147                raise ParsingError, msg
148           
149        # Fix some assertion failures
150        if(type(self.zone) == ArrayType and self.zone.shape == ()):
151            self.zone = self.zone[0]
152        if type(self.xllcorner) == ArrayType and self.xllcorner.shape == ():
153            self.xllcorner = self.xllcorner[0]
154        if type(self.yllcorner) == ArrayType and self.yllcorner.shape == ():
155            self.yllcorner = self.yllcorner[0]
156   
157        assert (type(self.xllcorner) == types.FloatType)
158        assert (type(self.yllcorner) == types.FloatType)
159        assert (type(self.zone) == types.IntType)
160       
161       
162    def change_points_geo_ref(self, points,
163                              points_geo_ref=None):
164        """
165        Change the geo reference of a list or Numeric array of points to
166        be this reference.(The reference used for this object)
167        If the points do not have a geo ref, assume 'absolute' values
168        """
169
170        is_list = False
171        if type(points) == types.ListType:
172            is_list = True
173
174        points = ensure_numeric(points, Float)
175       
176        if len(points.shape) == 1:
177            #One point has been passed
178            msg = 'Single point must have two elements'
179            assert len(points) == 2, msg
180            points = reshape(points, (1,2))
181
182        msg = 'Input must be an N x 2 array or list of (x,y) values. '
183        msg += 'I got an %d x %d array' %points.shape   
184        assert points.shape[1] == 2, msg               
185
186           
187        if points_geo_ref is not self:
188            #add point geo ref to points
189            if not points_geo_ref is None:
190                points[:,0] += points_geo_ref.xllcorner
191                points[:,1] += points_geo_ref.yllcorner
192       
193            #subtract primary geo ref from points
194            points[:,0] -= self.xllcorner
195            points[:,1] -= self.yllcorner
196
197           
198        if is_list:
199            points = points.tolist()
200           
201        return points
202
203   
204    def is_absolute(self):
205        """Return True if xllcorner==yllcorner==0 indicating that points
206        in question are absolute.
207        """
208
209        return allclose([self.xllcorner, self.yllcorner], 0) 
210
211       
212   
213    def get_absolute(self, points):
214        """
215        Given a set of points geo referenced to this instance,
216        return the points as absolute values.
217        """
218
219        #if self.is_absolute():
220        #    return points
221       
222
223        is_list = False
224        if type(points) == types.ListType:
225            is_list = True
226
227        points = ensure_numeric(points, Float)
228        if len(points.shape) == 1:
229            #One point has been passed
230            msg = 'Single point must have two elements'
231            if not len(points) == 2:
232                raise ShapeError, msg   
233                #points = reshape(points, (1,2))
234
235
236        msg = 'Input must be an N x 2 array or list of (x,y) values. '
237        msg += 'I got an %d x %d array' %points.shape   
238        if not points.shape[1] == 2:
239            raise ShapeError, msg   
240           
241       
242        # Add geo ref to points
243        if not self.is_absolute():
244            points[:,0] += self.xllcorner
245            points[:,1] += self.yllcorner
246
247       
248        if is_list:
249            points = points.tolist()
250             
251        return points
252
253
254    def reconcile_zones(self, other):
255        if other is None:
256            other = Geo_reference()
257        if self.zone == other.zone or\
258               self.zone == DEFAULT_ZONE and other.zone == DEFAULT_ZONE:
259            pass
260        elif self.zone == DEFAULT_ZONE:
261            self.zone = other.zone
262        elif other.zone == DEFAULT_ZONE:
263            other.zone = self.zone           
264        else:   
265            msg = 'Both geospatial_data objects must be in the same '+\
266                  'ZONE to allow reconciliation. I got zone %d and %d'\
267                  %(self.zone, other.zone)
268            raise ANUGAError, msg
269   
270    #def easting_northing2geo_reffed_point(self, x, y):
271    #    return [x-self.xllcorner, y - self.xllcorner]
272
273    #def easting_northing2geo_reffed_points(self, x, y):
274    #    return [x-self.xllcorner, y - self.xllcorner]
275
276    def get_origin(self):
277        return (self.zone, self.xllcorner, self.yllcorner)
278   
279    def __repr__(self):
280        return "(zone=%i easting=%f, northing=%f)" %(self.zone, self.xllcorner, self.yllcorner)
281   
282    def __cmp__(self,other):
283        # FIXME (DSG) add a tolerence
284        if other is None:
285            return 1
286        cmp = 0
287        if not (self.xllcorner == self.xllcorner):
288            cmp = 1
289        if not (self.yllcorner == self.yllcorner):
290            cmp = 1
291        if not (self.zone == self.zone):
292            cmp = 1
293        return cmp
294
295def write_NetCDF_georeference(origin, outfile):
296    """
297    Write georeferrence info to a netcdf file, usually sww.
298
299    The origin can be a georef instance or parrameters for a geo_ref instance
300
301    outfile is the name of the file to be written to.
302    """
303    if isinstance(origin, Geo_reference): 
304        geo_ref = origin
305    else:
306        geo_ref = apply(Geo_reference,origin)
307    geo_ref.write_NetCDF(outfile)
308    return geo_ref
309   
310#-----------------------------------------------------------------------
311
312if __name__ == "__main__":
313    pass
Note: See TracBrowser for help on using the repository browser.