source: inundation/coordinate_transforms/geo_reference.py @ 2959

Last change on this file since 2959 was 2942, checked in by duncan, 19 years ago

changing errors thrown by geo_ref. adding tests to geospatial
File size: 8.3 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
10from utilities.numerical_tools import ensure_numeric
11from utilities.anuga_exceptions import ANUGAError, TitleError, ParsingError
12
13
14#import exceptions
15#class TitleError(exceptions.IOError): pass
16#class ParsingError(exceptions.IOError): pass
17
18DEFAULT_ZONE = -1
19TITLE = '#geo reference' + "\n" #this title is referred to in the .xya format
20
21class Geo_reference:
22    """
23    """
24
25    def __init__(self,
26                 zone = DEFAULT_ZONE,
27                 xllcorner = 0.0,
28                 yllcorner = 0.0,
29                 datum = 'wgs84',
30                 projection = 'UTM',                 units = 'm',
31                 false_easting = 500000,
32                 false_northing = 10000000, #Default for southern hemisphere
33                 NetCDFObject=None,
34                 ASCIIFile=None,
35                 read_title=None):
36        """
37        input:
38        NetCDFObject - a handle to the netCDF file to be written to
39        ASCIIFile - a handle to the text file
40        read_title - the title of the georeference text, if it was read in.
41         If the function that calls this has already read the title line,
42         it can't unread it, so this info has to be passed.
43         If you know of a way to unread this info, then tell us.
44        """
45
46        self.false_easting = false_easting
47        self.false_northing = false_northing       
48        self.datum = datum
49        self.projection = projection
50        self.zone = zone       
51        self.units = units
52        self.xllcorner = xllcorner
53        self.yllcorner = yllcorner       
54           
55        if NetCDFObject is not None:
56            self.read_NetCDF(NetCDFObject)
57 
58        if ASCIIFile is not None:
59            self.read_ASCII(ASCIIFile, read_title=read_title)
60
61    def get_xllcorner(self):
62        return self.xllcorner
63   
64    def get_yllcorner(self):
65        return self.yllcorner
66   
67    def get_zone(self):
68        return self.zone
69       
70    def write_NetCDF(self, outfile):
71        outfile.xllcorner = self.xllcorner
72        outfile.yllcorner = self.yllcorner
73        outfile.zone = self.zone
74
75        outfile.false_easting = self.false_easting
76        outfile.false_northing = self.false_northing
77
78        outfile.datum = self.datum       
79        outfile.projection = self.projection
80        outfile.units = self.units
81
82    def read_NetCDF(self, infile):
83        self.xllcorner = infile.xllcorner[0]
84        self.yllcorner = infile.yllcorner[0] 
85        self.zone = infile.zone[0]
86
87        # Fix some assertion failures
88        if type(self.zone) == ArrayType and self.zone.shape == ():
89            self.zone = self.zone[0]
90        if type(self.xllcorner) == ArrayType and self.xllcorner.shape == ():
91            self.xllcorner = self.xllcorner[0]
92        if type(self.yllcorner) == ArrayType and self.yllcorner.shape == ():
93            self.yllcorner = self.yllcorner[0]
94
95        assert (type(self.xllcorner) == types.FloatType or\
96                type(self.xllcorner) == types.IntType)
97        assert (type(self.yllcorner) == types.FloatType or\
98                type(self.yllcorner) == types.IntType)
99        assert (type(self.zone) == types.IntType)
100
101        #FIXME (Ole) Read in the rest...
102       
103       
104    def write_ASCII(self, fd):
105        fd.write(TITLE)
106        fd.write(str(self.zone) + "\n") 
107        fd.write(str(self.xllcorner) + "\n") 
108        fd.write(str(self.yllcorner) + "\n")
109
110    def read_ASCII(self, fd,read_title=None):
111        try:
112            if read_title == None:
113                read_title = fd.readline() # remove the title line
114            if read_title[0:2].upper() != TITLE[0:2].upper():
115                msg = 'File error.  Expecting line: %s.  Got this line: %s' \
116                      %(TITLE, read_title)
117                raise TitleError, msg
118            self.zone = int(fd.readline())
119            self.xllcorner = float(fd.readline())
120            self.yllcorner = float(fd.readline())
121        except SyntaxError:
122                msg = 'File error.  Got syntax error while parsing geo reference' 
123                raise ParsingError, msg
124           
125        # Fix some assertion failures
126        if(type(self.zone) == ArrayType and self.zone.shape == ()):
127            self.zone = self.zone[0]
128        if type(self.xllcorner) == ArrayType and self.xllcorner.shape == ():
129            self.xllcorner = self.xllcorner[0]
130        if type(self.yllcorner) == ArrayType and self.yllcorner.shape == ():
131            self.yllcorner = self.yllcorner[0]
132   
133        assert (type(self.xllcorner) == types.FloatType)
134        assert (type(self.yllcorner) == types.FloatType)
135        assert (type(self.zone) == types.IntType)
136       
137        #false_easting = infile.false_easting[0]
138        #false_northing = infile.false_northing[0]
139       
140    def change_points_geo_ref(self, points,
141                              points_geo_ref=None):
142        """
143        Change the geo reference of a list or Numeric array of points to
144        be this reference.(The reference used for this object)
145        If the points do not have a geo ref, assume 'absolute' values
146        """
147
148        # FIXME: Use ensure_numeric
149       
150        is_list = False
151        if type(points) == types.ListType:
152            is_list = True
153
154        points = ensure_numeric(points, Float)
155       
156        if len(points.shape) == 1:
157            #One point has been passed
158            msg = 'Single point must have two elements'
159            assert len(points) == 2, msg
160            points = reshape(points, (1,2))
161
162           
163        if points_geo_ref is not self:
164            #add point geo ref to points
165            if not points_geo_ref is None:
166                points[:,0] += points_geo_ref.xllcorner
167                points[:,1] += points_geo_ref.yllcorner
168       
169            #subtract primary geo ref from points
170            points[:,0] -= self.xllcorner
171            points[:,1] -= self.yllcorner
172
173           
174        if is_list:
175            points = points.tolist()
176           
177        return points
178   
179   
180    def get_absolute(self, points):
181        """
182        Given a set of points geo referenced to this instance,
183        return the points as absolute values.
184        """
185       
186        is_list = False
187        if type(points) == types.ListType:
188            is_list = True
189            if len(points)>0 and type(points[0]) \
190                   not in [types.ListType,types.TupleType]:
191                #a single point is being passed.  make it a list of lists
192                points = [points]
193        elif type(points) == ArrayType:
194            if len(points.shape) == 1:
195                points = [points]       
196           
197        # convert into array
198        points = array(points).astype(Float) 
199        # add primary geo ref from points
200        points[:,0] += self.xllcorner
201        points[:,1] += self.yllcorner
202        if is_list:
203            points = points.tolist()
204             
205        return points
206
207
208    def reconcile_zones(self, other):
209
210        if self.zone == other.zone or self.zone == DEFAULT_ZONE and other.zone == DEFAULT_ZONE:
211            pass
212        elif self.zone == DEFAULT_ZONE:
213            self.zone = other.zone
214        elif other.zone == DEFAULT_ZONE:
215            other.zone = self.zone           
216        else:   
217            msg = 'Both geospatial_data objects must be in the same \
218            ZONE to allow reconciliation. I got zone %d and %d' %(self.zone, other.zone)
219            raise ANUGAError, msg
220   
221    #def easting_northing2geo_reffed_point(self, x, y):
222    #    return [x-self.xllcorner, y - self.xllcorner]
223
224    #def easting_northing2geo_reffed_points(self, x, y):
225    #    return [x-self.xllcorner, y - self.xllcorner]
226
227    def get_origin(self):
228        return (self.zone, self.xllcorner, self.yllcorner)
229   
230    def __repr__(self):
231        return "(zone=%i easting=%f, northing=%f)" %(self.zone, self.xllcorner, self.yllcorner)
232   
233    def __cmp__(self,other):
234        # FIXME (DSG) add a tolerence
235        if other is None:
236            return 1
237        cmp = 0
238        if not (self.xllcorner == self.xllcorner):
239            cmp = 1
240        if not (self.yllcorner == self.yllcorner):
241            cmp = 1
242        if not (self.zone == self.zone):
243            cmp = 1
244        return cmp
245       
246#-----------------------------------------------------------------------
247
248if __name__ == "__main__":
249    pass
Note: See TracBrowser for help on using the repository browser.