source: branches/inundation-numpy-branch/coordinate_transforms/geo_reference.py @ 6689

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

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