source: inundation/coordinate_transforms/geo_reference.py @ 2378

Last change on this file since 2378 was 2303, checked in by ole, 19 years ago

First stab at geo_spatial dataset and some tests
File size: 6.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
9from Numeric import array,Float,ArrayType
10
11DEFAULT_ZONE = 56
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       
78        assert (type(self.xllcorner) == types.FloatType or\
79                type(self.xllcorner) == types.IntType)
80        assert (type(self.yllcorner) == types.FloatType or\
81                type(self.yllcorner) == types.IntType)
82        assert (type(self.zone) == types.IntType)
83
84        #FIXME (Ole) Read in the rest...
85       
86       
87    def write_ASCII(self, fd):
88        fd.write(TITLE)
89        fd.write(str(self.zone) + "\n") 
90        fd.write(str(self.xllcorner) + "\n") 
91        fd.write(str(self.yllcorner) + "\n")
92
93    def read_ASCII(self, fd,read_title=None):
94        if read_title == None:
95            read_title = fd.readline() # remove the title line
96        #print "gr read_title[0:2].upper()",read_title[0:2].upper()
97        #print "gr TITLE[0:2].upper()",TITLE[0:2].upper()
98        if read_title[0:2].upper() != TITLE[0:2].upper():
99            raise SyntaxError
100        self.zone = int(fd.readline())
101        self.xllcorner = float(fd.readline())
102        self.yllcorner = float(fd.readline())
103           
104       
105        assert (type(self.xllcorner) == types.FloatType)
106        assert (type(self.yllcorner) == types.FloatType)
107        assert (type(self.zone) == types.IntType)
108       
109        #false_easting = infile.false_easting[0]
110        #false_northing = infile.false_northing[0]
111       
112    def change_points_geo_ref(self, points,
113                              points_geo_ref=None):
114        """
115        Change the geo reference of a list or Numeric array of points to
116        be this reference.(The reference used for this object)
117        If the points do not have a geo ref, assume 'absolute' values
118        """
119       
120        is_list = False
121        if type(points) == types.ListType:
122            is_list = True
123            if len(points)>0 and type(points[0]) not in [types.ListType,types.TupleType]:
124                #a single point is being passed.  make it a list of lists
125                points = [points]
126        elif type(points) == ArrayType:
127            if len(points.shape) == 1:
128                points = [points]
129       
130        # currently untested
131        # Always returns a list of lists
132        if points_geo_ref is self:
133            return points       
134           
135        #convert into array
136        points = array(points).astype(Float) 
137        #add point geo ref to points
138        if not points_geo_ref is None:
139            points[:,0] += points_geo_ref.xllcorner
140            points[:,1] += points_geo_ref.yllcorner
141       
142        #subtract primary geo ref from points
143        points[:,0] -= self.xllcorner
144        points[:,1] -= self.yllcorner
145        if is_list:
146            points = points.tolist()
147        return points
148   
149    def get_absolute(self, points):
150        """
151        Given a set of points geo referenced to this instance,
152        return the points as absolute values.
153        """
154       
155        is_list = False
156        if type(points) == types.ListType:
157            is_list = True
158            if len(points)>0 and type(points[0]) \
159                   not in [types.ListType,types.TupleType]:
160                #a single point is being passed.  make it a list of lists
161                points = [points]
162        elif type(points) == ArrayType:
163            if len(points.shape) == 1:
164                points = [points]       
165           
166        # convert into array
167        points = array(points).astype(Float) 
168        # add primary geo ref from points
169        points[:,0] += self.xllcorner
170        points[:,1] += self.yllcorner
171        if is_list:
172            points = points.tolist()
173             
174        return points
175
176   
177    #def easting_northing2geo_reffed_point(self, x, y):
178    #    return [x-self.xllcorner, y - self.xllcorner]
179
180    #def easting_northing2geo_reffed_points(self, x, y):
181    #    return [x-self.xllcorner, y - self.xllcorner]
182
183    def get_origin(self):
184        return (self.zone, self.xllcorner, self.yllcorner)
185   
186    def __repr__(self):
187        return "(zone=%i easting=%f, northing=%f)" %(self.zone, self.xllcorner, self.yllcorner)
188   
189    def __cmp__(self,other):
190        # FIXME (DSG) add a tolerence
191        if other is None:
192            return 1
193        cmp = 0
194        if not (self.xllcorner == self.xllcorner):
195            cmp = 1
196        if not (self.yllcorner == self.yllcorner):
197            cmp = 1
198        if not (self.zone == self.zone):
199            cmp = 1
200        return cmp
201       
202#-----------------------------------------------------------------------
203
204if __name__ == "__main__":
205    pass
Note: See TracBrowser for help on using the repository browser.