source: inundation/coordinate_transforms/geo_reference.py @ 2608

Last change on this file since 2608 was 2608, checked in by ole, 18 years ago

Played with custom exceptions for ANUGA

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