/[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 4443 - (show annotations)
Mon Jun 10 04:50:27 2013 UTC (6 years, 4 months ago) by caltinay
File MIME type: text/x-python
File size: 14282 byte(s)
Update test to take into account changed height scale.

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 import esys.escript as esc
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
71 .. note:: needs to be overwritten by a particular reference system
72 """
73 raise NotImplementedError()
74
75 def isCartesian(self):
76 """
77 returns if the reference system is Cartesian
78
79 .. note:: needs to be overwritten by a particular reference system
80
81 :rtype: ``bool``
82 """
83 raise NotImplementedError()
84
85 def createTransformation(self, domain):
86 """
87 creates an appropriate coordinate transformation on a given domain
88
89 .. note:: needs to be overwritten by a particular reference system
90
91 :param domain: domain of transformation
92 :type domain: `esys.escript.AbstractDomain`
93 :rtype: `SpatialCoordinateTransformation`
94 """
95 raise NotImplementedError()
96
97 class CartesianReferenceSystem(ReferenceSystem):
98 """
99 Identifies the Cartesian coordinate system
100 """
101 def __init__(self, name="CARTESIAN"):
102 """
103 set up Cartesian coordinate system
104 """
105 super(CartesianReferenceSystem, self).__init__(name)
106
107 def isTheSame(self, other):
108 """
109 test if argument ``other`` defines the same reference system
110
111 :param other: a second reference system
112 :type other: `ReferenceSystem`
113 :returns: ``True`` if ``other`` is a `CartesianReferenceSystem` instance.
114 :rtype: ``bool``
115 :note: every two `CartesianReferenceSystem` instances are considered
116 as being the same.
117 """
118 return isinstance(other, CartesianReferenceSystem)
119
120 def createTransformation(self, domain):
121 """
122 creates an appropriate coordinate transformation on a given domain
123
124 :param domain: domain of transformation
125 :type domain: `esys.escript.AbstractDomain`
126 :rtype: `SpatialCoordinateTransformation`
127 """
128 return SpatialCoordinateTransformation(domain, reference=self)
129
130 def isCartesian(self):
131 """
132 returns if the reference system is Cartesian
133
134 :rtype: ``bool``
135 """
136 return True
137
138 class GeodeticReferenceSystem(ReferenceSystem):
139 """
140 Identifies a Geodetic coordinate system
141 """
142 def __init__(self, a=6378137.0 *U.m, f=1/298.257223563, angular_unit=1.*U.DEG, height_unit=1.*U.km, name="WGS84"):
143 """
144 initializes a geodetic reference system
145
146 :param a: semi-major axis in meter
147 :type a: positive ``double``
148 :param f: flattening
149 :type f: non-negative ``double``, less than one
150 :param name: name of the reference system
151 :type name: ``str``
152 :param angular_unit: factor to scale the unit of latitude and
153 longitude to radians.
154 :type angular_unit: positive ``double``
155 :param height_unit: factor to scale the unit of latitude and
156 longitude to radians.
157 :type height_unit: positive ``double``
158 """
159 if not a>0:
160 raise ValueError("length of semi-major axis a must be positive.")
161 if not ( f>=0 and f<1 ):
162 raise ValueError("flattening f must be non-negative and less than one.")
163 if not angular_unit > 0:
164 raise ValueError("angular_unit must be positive.")
165 if not height_unit > 0:
166 raise ValueError("height_unit must be positive.")
167 super(GeodeticReferenceSystem, self).__init__(name)
168
169 self.__a=a
170 self.__f=f
171 self.__angular_unit=angular_unit
172 self.__height_unit=height_unit
173
174
175 def isCartesian(self):
176 """
177 returns if the reference system is Cartesian
178
179 :rtype: ``bool``
180 """
181 return False
182
183 def getAngularUnit(self):
184 """
185 returns the angular unit
186 """
187 return self.__angular_unit
188
189 def getHeightUnit(self):
190 """
191 returns the height unit
192 """
193 return self.__height_unit
194
195 def getSemiMajorAxis(self):
196 """
197 returns the length of semi major axis
198 """
199 return self.__a
200
201 def getSemiMinorAxis(self):
202 """
203 returns the length of semi minor axis
204 """
205 a=self.getSemiMajorAxis()
206 f=self.getFlattening()
207 return a*(1-f)
208
209 def getFlattening(self):
210 """
211 returns the flattening
212 """
213 return self.__f
214
215 def isTheSame(self, other):
216 """
217 test if ``other`` defines the same reference system
218
219 :param other: a second reference system
220 :type other: `ReferenceSystem`
221 :returns: ``True`` if other defines then same reference system
222 :rtype: ``bool``
223
224 .. note:: two `GeodeticReferenceSystem` are considered to be the same
225 if the use the same semi major axis, the same flattening
226 and the same angular unit.
227 """
228 if isinstance(other, GeodeticReferenceSystem):
229 if self.getSemiMajorAxis() == other.getSemiMajorAxis() \
230 and self.getFlattening() == other.getFlattening() \
231 and self.getAngularUnit() == other.getAngularUnit():
232 return True
233 else:
234 return False
235 else:
236 return False
237
238 def createTransformation(self, domain):
239 """
240 creates an appropriate coordinate transformation on a given domain
241
242 :param domain: domain of transformation
243 :type domain: `esys.escript.AbstractDomain`
244 :rtype: `SpatialCoordinateTransformation`
245 """
246 return GeodeticCoordinateTransformation(domain, reference=self)
247
248 def SphericalReferenceSystem(R=6378137.0*U.m):
249 """
250 returns the `GeodeticReferenceSystem` of a sphere
251 :param R: sphere radius
252 :type R: positive ``double``
253 """
254 return GeodeticReferenceSystem(a=R, f=0, angular_unit=1*U.DEG, height_unit=1.*U.km, name="SPHERE")
255
256 def WGS84ReferenceSystem():
257 """
258 returns the `GeodeticReferenceSystem` for the WGS84 Ellipsoid
259 """
260 return GeodeticReferenceSystem(a=6378137.0 *U.m, f=1/298.257223563, angular_unit=1*U.DEG, height_unit=100.*U.km, name="WGS84")
261
262 def GRS80ReferenceSystem():
263 """
264 returns the `GeodeticReferenceSystem` for the GRS80 Ellipsoid eg. used by Geocentric Datum of Australia GDA94
265 """
266 return GeodeticReferenceSystem(a=6378137.0 *U.m, f=1/298.257222101, angular_unit=1*U.DEG, height_unit=1.*U.km, name="GRS80")
267
268
269 class SpatialCoordinateTransformation(object):
270 """
271 Defines an orthogonal coordinate transformation from a domain into the
272 Cartesian domain using a coordinate transformation.
273
274 The default implementation is the identity transformation (i.e.
275 no changes are applied to the domain). Overwrite the appropriate
276 methods to define other coordinate systems.
277 """
278 def __init__(self, domain, reference=CartesianReferenceSystem()):
279 """
280 set up the orthogonal coordinate transformation.
281
282 :param domain: domain in the domain of the coordinate transformation
283 :type domain: `esys.escript.AbstractDomain`
284 :param reference: the reference system
285 :type reference: `ReferenceSystem`
286 """
287 self.__domain = domain
288 self.__reference_system=reference
289 self._volumefactor=esc.Scalar(1., esc.Function(domain))
290 self._scaling_factors = esc.Vector(1., esc.Function(domain))
291
292 def __eq__(self, other):
293 return self.isTheSame(other)
294
295 def __ne__(self, other):
296 return not self.isTheSame(other)
297
298 def isTheSame(self, other):
299 """
300 test if argument ``other`` defines the same coordinate transformation
301
302 :param other: a second coordinate transformation
303 :type other: `SpatialCoordinateTransformation`
304 :returns: ``True`` if other defines then same coordinate transformation
305 :rtype: ``bool``
306 """
307 if isinstance(other, SpatialCoordinateTransformation):
308 if self.getDomain() == other.getDomain() \
309 and self.getReferenceSystem() == other.getReferenceSystem():
310 return True
311
312 return False
313
314 def isCartesian(self):
315 """
316 returns ``True`` if the scaling factors (and the volume factor) are equal to 1
317
318 :rtype: ``bool``
319 """
320 return self.__reference_system.isCartesian()
321
322 def getDomain(self):
323 """
324 returns the domain of the coordinate transformation.
325
326 :rtype: `esys.escript.AbstractDomain`
327 """
328 return self.__domain
329
330 def getReferenceSystem(self):
331 """
332 returns the reference system used to to define the coordinate transformation
333
334 :rtype: `ReferenceSystem`
335 """
336 return self.__reference_system
337
338 def getVolumeFactor(self):
339 """
340 returns the volume factor for the coordinate transformation
341
342 :rtype: `esys.escript.Scalar`
343 """
344 return self._volumefactor
345
346 def getScalingFactors(self):
347 """
348 returns the scaling factors
349
350 :rtype: `esys.escript.Vector`
351 """
352 return self._scaling_factors
353
354 def CartesianCoordinateTransformation(domain, reference=CartesianReferenceSystem() ):
355
356 return SpatialCoordinateTransformation(domain, reference)
357
358 class GeodeticCoordinateTransformation(SpatialCoordinateTransformation):
359 """
360 A geodetic coordinate transformation
361 """
362 def __init__(self, domain, reference=WGS84ReferenceSystem() ):
363 """
364 set up the orthogonal coordinate transformation.
365
366 :param domain: domain in the domain of the coordinate transformation
367 :type domain: `esys.escript.AbstractDomain`
368 :param reference: the reference system
369 :type reference: `ReferenceSystem`
370
371 """
372 DIM=domain.getDim()
373 super(GeodeticCoordinateTransformation, self).__init__(domain, reference )
374
375 a=reference.getSemiMajorAxis()
376 f=reference.getFlattening()
377 f_a=reference.getAngularUnit()
378 f_h=reference.getHeightUnit()
379
380 x=esc.Function(domain).getX()
381 if DIM == 2:
382 phi=0.
383 else:
384 phi=x[1] * f_a
385 h=x[DIM-1] * f_h
386
387 e = esc.sqrt(2*f-f**2)
388 N = a/esc.sqrt(1 - e**2 * esc.sin(phi)**2 )
389 M = ( a*(1-e**2) ) /esc.sqrt(1 - e**2 * esc.sin(phi)**2 )**3
390 v_phi = f_a * (M + h)
391 v_lam = f_a * (N + h) * esc.cos(phi)
392 v_h = f_h
393 s= esc.Vector(1., esc.Function(domain))
394 if DIM == 2:
395 v= v_phi * v_h
396 s[0]=1/v_lam
397 s[1]=1/v_h
398 else:
399 v= v_phi * v_lam * v_h
400 s[0]=1/v_lam
401 s[1]=1/v_phi
402 s[2]=1/v_h
403
404 self._volumefactor=v
405 self._scaling_factors = s
406
407 def makeTranformation(domain, coordinates=None):
408 """
409 returns a `SpatialCoordinateTransformation` for the given domain
410
411 :param domain: domain in the domain of the coordinate transformation
412 :type domain: `esys.escript.AbstractDomain`
413 :param coordinates: the reference system or spatial coordinate system.
414 :type coordinates: `ReferenceSystem` or `SpatialCoordinateTransformation`
415 :return: the spatial coordinate system for the given domain of the specified
416 reference system ``coordinates``. If ``coordinates`` is already spatial coordinate system based on the
417 riven domain ``coordinates`` is returned. Otherwise an appropriate spatial coordinate system
418 is created.
419 :rtype: `SpatialCoordinateTransformation`
420 """
421 if coordinates == None:
422 return CartesianCoordinateTransformation(domain)
423 elif isinstance(coordinates, ReferenceSystem):
424 return coordinates.createTransformation(domain)
425 else:
426 if not coordinates.getDomain() == domain:
427 raise ValueError("Domain of spatial coordinate system and given domain don't match.")
428 else:
429 return coordinates
430

  ViewVC Help
Powered by ViewVC 1.1.26