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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4433 - (hide annotations)
Fri May 31 12:09:58 2013 UTC (5 years, 11 months ago) by gross
File MIME type: text/x-python
File size: 14266 byte(s)
some clarifications on geodetic coordinates. 
order of background magnetic flux density component has been corrected: input is now B_east, B_north, B_vertical.


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

  ViewVC Help
Powered by ViewVC 1.1.26