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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 4432 by jfenwick, Thu May 9 08:42:44 2013 UTC revision 4433 by gross, Fri May 31 12:09:58 2013 UTC
# Line 139  class GeodeticReferenceSystem(ReferenceS Line 139  class GeodeticReferenceSystem(ReferenceS
139      """      """
140      Identifies a Geodetic coordinate system      Identifies a Geodetic coordinate system
141      """      """
142      def __init__(self, a=6378137.0 *U.m, f=1/298.257223563, angular_unit=1.*U.DEG, name="WGS84"):      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          initializes a geodetic reference system
145    
# Line 152  class GeodeticReferenceSystem(ReferenceS Line 152  class GeodeticReferenceSystem(ReferenceS
152          :param angular_unit: factor to scale the unit of latitude and          :param angular_unit: factor to scale the unit of latitude and
153                               longitude to radians.                               longitude to radians.
154          :type  angular_unit: positive ``double``          :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:          if not a>0:
160              raise ValueError("length of semi-major axis a must be positive.")              raise ValueError("length of semi-major axis a must be positive.")
# Line 159  class GeodeticReferenceSystem(ReferenceS Line 162  class GeodeticReferenceSystem(ReferenceS
162              raise ValueError("flattening f must be non-negative and less than one.")              raise ValueError("flattening f must be non-negative and less than one.")
163          if not angular_unit > 0:          if not angular_unit > 0:
164              raise ValueError("angular_unit must be positive.")              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)          super(GeodeticReferenceSystem, self).__init__(name)
168          
169          self.__a=a          self.__a=a
170          self.__f=f          self.__f=f
171          self.__angular_unit=angular_unit          self.__angular_unit=angular_unit
172            self.__height_unit=height_unit
173            
174    
175      def isCartesian(self):      def isCartesian(self):
176          """          """
# Line 177  class GeodeticReferenceSystem(ReferenceS Line 185  class GeodeticReferenceSystem(ReferenceS
185          returns the angular unit          returns the angular unit
186          """          """
187          return self.__angular_unit          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):      def getSemiMajorAxis(self):
196          """          """
197          returns the length of semi major axis          returns the length of semi major axis
# Line 237  def SphericalReferenceSystem(R=6378137.0 Line 251  def SphericalReferenceSystem(R=6378137.0
251      :param R: sphere radius      :param R: sphere radius
252      :type R: positive ``double``      :type R: positive ``double``
253      """      """
254      return GeodeticReferenceSystem(a=R, f=0, angular_unit=1*U.DEG, name="SPHERE")      return GeodeticReferenceSystem(a=R, f=0, angular_unit=1*U.DEG,  height_unit=1.*U.km, name="SPHERE")
255    
256  def WGS84ReferenceSystem():  def WGS84ReferenceSystem():
257      """      """
258      returns the `GeodeticReferenceSystem` for the WGS84 Ellipsoid      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, name="WGS84")      return GeodeticReferenceSystem(a=6378137.0 *U.m, f=1/298.257223563, angular_unit=1*U.DEG,  height_unit=1.*U.km, name="WGS84")
261    
262  def GRS80ReferenceSystem():  def GRS80ReferenceSystem():
263      """      """
264      returns the `GeodeticReferenceSystem` for the GRS80 Ellipsoid eg. used by Geocentric Datum of Australia GDA94      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, name="GRS80")      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):  class SpatialCoordinateTransformation(object):
# Line 360  class GeodeticCoordinateTransformation(S Line 374  class GeodeticCoordinateTransformation(S
374    
375          a=reference.getSemiMajorAxis()          a=reference.getSemiMajorAxis()
376          f=reference.getFlattening()          f=reference.getFlattening()
377            f_a=reference.getAngularUnit()
378            f_h=reference.getHeightUnit()
379    
380          x=esc.Function(domain).getX()          x=esc.Function(domain).getX()
381          phi=x[0] * reference.getAngularUnit()          if DIM == 2:
382          h=x[DIM-1]         phi=0.
383        else:
384               phi=x[1] * f_a
385            h=x[DIM-1] * f_h
386    
387            
388          e = esc.sqrt(2*f-f**2)          e = esc.sqrt(2*f-f**2)
389          N = a/esc.sqrt(1 - e**2 * esc.sin(phi)**2 )          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          M = ( a*(1-e**2) ) /esc.sqrt(1 - e**2 * esc.sin(phi)**2 )**3
391          v_phi = (M + h)          v_phi = f_a * (M + h)
392          v_lam = (N + h) * esc.cos(phi)          v_lam = f_a * (N + h) * esc.cos(phi)
393            v_h = f_h
394          s= esc.Vector(1., esc.Function(domain))          s= esc.Vector(1., esc.Function(domain))
395          if DIM == 2:          if DIM == 2:
396              v= v_phi              v= v_phi * v_h
397              s=esc.Vector(1., esc.Function(domain))              s[0]=1/v_lam
398              s[0]=1/v_phi              s[1]=1/v_h
399          else:          else:
400              v= v_phi * v_lam              v= v_phi * v_lam * v_h
401              s[0]=1/v_phi              s[0]=1/v_lam
402              s[1]=1/v_lam              s[1]=1/v_phi
403                s[2]=1/v_h
404    
405          self._volumefactor=v          self._volumefactor=v
406          self._scaling_factors = s          self._scaling_factors = s
407    

Legend:
Removed from v.4432  
changed lines
  Added in v.4433

  ViewVC Help
Powered by ViewVC 1.1.26