/[escript]/trunk/downunder/py_src/coordinates.py
ViewVC logotype

Contents of /trunk/downunder/py_src/coordinates.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4389 - (show annotations)
Mon May 6 01:52:05 2013 UTC (6 years, 5 months ago) by caltinay
File MIME type: text/x-python
File size: 13392 byte(s)
retab coordinates.py

1
2 ##############################################################################
3 #
4 # Copyright (c) 2003-2013 by University of Queensland
5 # http://www.uq.edu.au
6 #
7 # Primary Business: Queensland, Australia
8 # Licensed under the Open Software License version 3.0
9 # http://www.opensource.org/licenses/osl-3.0.php
10 #
11 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 # Development since 2012 by School of Earth Sciences
13 #
14 ##############################################################################
15
16 """Functions to deal with coordinate systems"""
17
18 __copyright__="""Copyright (c) 2003-2013 by University of Queensland
19 http://www.uq.edu.au
20 Primary Business: Queensland, Australia"""
21 __license__="""Licensed under the Open Software License version 3.0
22 http://www.opensource.org/licenses/osl-3.0.php"""
23 __url__="https://launchpad.net/escript-finley"
24
25 __all__ = ['ReferenceSystem', 'CartesianReferenceSystem',
26 'GeodeticReferenceSystem', 'SphericalReferenceSystem',
27 'WGS84ReferenceSystem', 'GRS80ReferenceSystem',
28 'SpatialCoordinateTransformation', 'GeodeticCoordinateTransformation',
29 'CartesianCoordinateTransformation', 'makeTranformation']
30
31 from esys.escript import unitsSI as U
32 from esys.escript import *
33
34 class ReferenceSystem(object):
35 """
36 Generic identifier for coordinate systems.
37 """
38 def __init__(self,name='none'):
39 """
40 initialization of reference system
41
42 :param name: name of the reference system
43 :type name : ``str``
44 """
45 self.__name=name
46
47 def getName(self):
48 """
49 returns the name of the reference system
50 """
51 return self.__name
52
53 def __str__(self):
54 return ("%s (id %s)"%(self.getName(),id(self)))
55
56 def __eq__(self, other):
57 return self.isTheSame(other)
58
59 def __ne__(self, other):
60 return not self.isTheSame(other)
61
62 def isTheSame(self, other):
63 """
64 test if argument ``other`` defines the same reference system
65
66 :param other: a second reference system
67 :type other: `ReferenceSystem`
68 :returns: ``True`` if other defines the same reference system
69 :rtype: ``bool``
70 :note: needs to be overwritten by a particular reference system
71 """
72 raise NotImplementedError()
73
74 def isCartesian(self):
75 """
76 returns if the reference system is Cartesian
77 :rtype: ``bool``
78 \ :note: needs to be overwritten by a particular reference system
79 """
80 raise NotImplementedError()
81
82 def createTransformation(self, domain):
83 """
84 creates an appropriate coordinate transformation on a given domain
85
86 :param domain: domain of transformation
87 :type domain: `esys.escript.AbstractDomain'
88 :rtype: `SpatialCoordinateTransformation`
89 :note: needs to be overwritten by a particular reference system
90 """
91 raise NotImplementedError()
92
93 class CartesianReferenceSystem(ReferenceSystem):
94 """
95 Identifies the Cartesian coordinate system
96 """
97 def __init__(self, name="CARTESIAN"):
98 """
99 set up Cartesian coordinate system
100 """
101 super(CartesianReferenceSystem, self).__init__(name)
102
103 def isTheSame(self, other):
104 """
105 test if argument ``other`` defines the same reference system
106
107 :param other: a second reference system
108 :type other: `ReferenceSystem`
109 :returns: ``True`` if ``other`` is a `CartesianReferenceSystem` instance.
110 :rtype: ``bool``
111 :note: every two `CartesianReferenceSystem` instances are considered
112 as being the same.
113 """
114 return isinstance(other, CartesianReferenceSystem)
115
116 def createTransformation(self, domain):
117 """
118 creates an appropriate coordinate transformation on a given domain
119
120 :param domain: domain of transformation
121 :type domain: `esys.escript.AbstractDomain'
122 :rtype: `SpatialCoordinateTransformation`
123 """
124 return SpatialCoordinateTransformation(domain, reference=self)
125
126 def isCartesian(self):
127 """
128 returns if the reference system is Cartesian
129 :rtype: ``bool``
130 """
131 return True
132
133 class GeodeticReferenceSystem(ReferenceSystem):
134 """
135 Identifies a Geodetic coordinate system
136 """
137 def __init__(self, a=6378137.0 *U.m, f=1/298.257223563, angular_unit=1.*U.DEG, name="WGS84"):
138 """
139 initializes a geodetic reference system
140
141 :param a: semi-major axis in meter
142 :type a: positive ``double``
143 :param f: flattening
144 :type f: non-negative ``double``, less than one
145 :param name: name of the reference system
146 :type name: ``str``
147 :param angular_unit: factor to scale the unit of latitude and
148 longitude to radians.
149 :type angular_unit: positive ``double``
150 """
151 if not a>0:
152 raise ValueError("length of semi-major axis a must be positive.")
153 if not ( f>=0 and f<1 ):
154 raise ValueError("flattening f must be non-negative and less than one.")
155 if not angular_unit > 0:
156 raise ValueError("angular_unit must be positive.")
157 super(GeodeticReferenceSystem, self).__init__(name)
158 self.__a=a
159 self.__f=f
160 self.__angular_unit=angular_unit
161
162 def isCartesian(self):
163 """
164 returns if the reference system is Cartesian
165 :rtype: ``bool``
166 """
167 return False
168
169 def getAngularUnit(self):
170 """
171 returns the angular unit
172 """
173 return self.__angular_unit
174
175 def getSemiMajorAxis(self):
176 """
177 returns the length of semi major axis
178 """
179 return self.__a
180
181 def getSemiMinorAxis(self):
182 """
183 returns the length of semi minor axis
184 """
185 a=self.getSemiMajorAxis()
186 f=self.getFlattening()
187 return a*(1-f)
188
189 def getFlattening(self):
190 """
191 returns the flattening
192 """
193 return self.__f
194
195 def isTheSame(self, other):
196 """
197 test if ``other`` defines the same reference system
198
199 :param other: a second reference system
200 :type other: `ReferenceSystem`
201 :returns: ``True`` if other defines then same reference system
202 :rtype: ``bool``
203 :note: two `GeodeticReferenceSystem` are considered to be the same
204 if the use the same semi major axis, the same flattening
205 and the same angular unit.
206 """
207 if isinstance(other, GeodeticReferenceSystem):
208 if self.getSemiMajorAxis() == other.getSemiMajorAxis() \
209 and self.getFlattening() == other.getFlattening() \
210 and self.getAngularUnit() == other.getAngularUnit():
211 return True
212 else:
213 return False
214 else:
215 return False
216
217 def createTransformation(self, domain):
218 """
219 creates an appropriate coordinate transformation on a given domain
220
221 :param domain: domain of transformation
222 :type domain:`esys.escript.AbstractDomain'
223 :rtype: `SpatialCoordinateTransformation`
224 """
225 return GeodeticCoordinateTransformation(domain, reference=self)
226
227 def SphericalReferenceSystem(R=6378137.0*U.m):
228 """
229 returns the `GeodeticReferenceSystem` of a sphere
230 :param R: sphere radius
231 :type R: positive ``double``
232 """
233 return GeodeticReferenceSystem(a=R, f=0, angular_unit=1*U.DEG, name="SPHERE")
234
235 def WGS84ReferenceSystem():
236 """
237 returns the `GeodeticReferenceSystem` for the WGS84 Ellipsoid
238 """
239 return GeodeticReferenceSystem(a=6378137.0 *U.m, f=1/298.257223563, angular_unit=1*U.DEG, name="WGS84")
240
241 def GRS80ReferenceSystem():
242 """
243 returns the `GeodeticReferenceSystem` for the GRS80 Ellipsoid eg. used by Geocentric Datum of Australia GDA94
244 """
245 return GeodeticReferenceSystem(a=6378137.0 *U.m, f=1/298.257222101, angular_unit=1*U.DEG, name="GRS80")
246
247
248 class SpatialCoordinateTransformation(object):
249 """
250 Defines an orthogonal coordinate transformation from a domain into the
251 Cartesian domain using a coordinate transformation.
252
253 The default implementation is the identity transformation (i.e.
254 no changes are applied to the domain). Overwrite the appropriate
255 methods to define other coordinate systems.
256 """
257 def __init__(self, domain, reference=CartesianReferenceSystem()):
258 """
259 set up the orthogonal coordinate transformation.
260
261 :param domain: domain in the domain of the coordinate transformation
262 :type domain: `esys.escript.AbstractDomain'
263 :param reference: the reference system
264 :type reference: `ReferenceSystem`
265 """
266 self.__domain = domain
267 self.__reference_system=reference
268 self._volumefactor=Scalar(1., Function(domain))
269 self._scaling_factors = Vector(1., Function(domain))
270
271 def __eq__(self, other):
272 return self.isTheSame(other)
273
274 def __ne__(self, other):
275 return not self.isTheSame(other)
276
277 def isTheSame(self, other):
278 """
279 test if argument ``other`` defines the same coordinate transformation
280
281 :param other: a second coordinate transformation
282 :type other: `SpatialCoordinateTransformation`
283 :returns: ``True`` if other defines then same coordinate transformation
284 :rtype: ``bool``
285 """
286 if isinstance(other, SpatialCoordinateTransformation):
287 if self.getDomain() == other.getDomain() \
288 and self.getReferenceSystem() == other.getReferenceSystem():
289 return True
290
291 return False
292
293 def isCartesian(self):
294 """
295 returns ``True`` if the scaling factors (and the volume factor) are equal to 1
296 :rtype: ``bool``
297 """
298 return self.__reference_system.isCartesian()
299
300 def getDomain(self):
301 """
302 returns the domain of the coordinate transformation.
303 :rtype: `esys.escript.AbstractDomain'
304 """
305 return self.__domain
306
307 def getReferenceSystem(self):
308 """
309 returns the reference system used to to define the coordinate transformation
310 :rtype: `ReferenceSystem`
311 """
312 return self.__reference_system
313
314 def getVolumeFactor(self):
315 """
316 returns the volume factor for the coordinate transformation
317 :rtype: `esys.escript.Scalar'
318 """
319 return self._volumefactor
320
321 def getScalingFactors(self):
322 """
323 returns the scaling factors
324 :rtype: `esys.escript.Vector'
325 """
326 return self._scaling_factors
327
328 def CartesianCoordinateTransformation(domain, reference=CartesianReferenceSystem() ):
329
330 return SpatialCoordinateTransformation(domain, reference)
331
332 class GeodeticCoordinateTransformation(SpatialCoordinateTransformation):
333 """
334 A geodetic coordinate transformation
335 """
336 def __init__(self, domain, reference=WGS84ReferenceSystem() ):
337 """
338 set up the orthogonal coordinate transformation.
339
340 :param domain: domain in the domain of the coordinate transformation
341 :type domain: `esys.escript.AbstractDomain'
342 :param reference: the reference system
343 :type reference: `ReferenceSystem`
344
345 """
346 DIM=domain.getDim()
347 super(GeodeticCoordinateTransformation, self).__init__(domain, reference )
348
349 a=reference.getSemiMajorAxis()
350 f=reference.getFlattening()
351
352 x=Function(domain).getX()
353 phi=x[0] * reference.getAngularUnit()
354 h=x[DIM-1]
355
356 e = sqrt(2*f-f**2)
357 N = a/sqrt(1 - e**2 * sin(phi)**2 )
358 M = ( a*(1-e**2) ) /sqrt(1 - e**2 * sin(phi)**2 )**3
359 v_phi = (M + h)
360 v_lam = (N + h) * cos(phi)
361 s= Vector(1., Function(domain))
362 if DIM == 2:
363 v= v_phi
364 s=Vector(1., Function(domain))
365 s[0]=1/v_phi
366 else:
367 v= v_phi * v_lam
368 s[0]=1/v_phi
369 s[1]=1/v_lam
370 self._volumefactor=v
371 self._scaling_factors = s
372
373 def makeTranformation(domain, coordinates=None):
374 """
375 returns a `SpatialCoordinateTransformation` for the given domain
376
377 :param domain: domain in the domain of the coordinate transformation
378 :type domain: `esys.escript.AbstractDomain'
379 :param coordinates: the reference system or spatial coordinate system.
380 :type coordinates: `ReferenceSystem` or `SpatialCoordinateTransformation`
381 :return: the spatial coordinate system for the given domain of the specified
382 reference system ``coordinates``. If ``coordinates`` is already spatial coordinate system based on the
383 riven domain ``coordinates`` is returned. Otherwise an appropriate spatial coordinate system
384 is created.
385 :rtype: `SpatialCoordinateTransformation`
386 """
387 if coordinates == None:
388 return CartesianCoordinateTransformation(domain)
389 elif isinstance(coordinates, ReferenceSystem):
390 return coordinates.createTransformation(domain)
391 else:
392 if not coordinates.getDomain() == domain:
393 raise ValueError("Domain of spatial coordinate system and given domain don't match.")
394 else:
395 return coordinates
396

  ViewVC Help
Powered by ViewVC 1.1.26