source: inundation/coordinate_transforms/geo_reference.py @ 2593

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

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