Così ho iniziato con il tutorial e con Dive Into Python.
Dopo qualche esercizio ho avuto bisogno di qualcosa di più concreto da fare: quale occasione migliore per ripescare nella memoria qualche confuso ricordo di cartografia?
Quel che segue è un brutto codice, sicuramente ridondante... benvenuti i suggerimenti e le correzioni!
Quasi tutte le fonti sono citate nei commenti.
#! /usr/bin/env python
# -*- coding: utf-8 -*-
from __future__ import division
import math
class ellipsoid():
"""
b semi-minor axis (polar)
a semi-major axis (equatorial)
f flattering
rf reverse (inverse) flattering
"""
def __init__(self, a,**kargs):
self.a=a
self.b=None
## a sloppy and shallow parsing...
if len(kargs)<1:
exit
elif 'b' in kargs.keys():
self.b = kargs['b']
elif 'f' in kargs.keys():
self.b =self.a*(1-kargs['f'])
elif 'rf' in kargs.keys():
self.b=self.a*(1-1/kargs['rf'])
if 'name' in kargs.keys():
self.name=kargs['name']
else:
self.name=None
def getFirstEccentricity(self):
"""
first eccentricity = sqrt(1-(b2/a2)**2)
"""
return math.sqrt(1 - (self.b / self.a) ** 2)
def getFirstEccentricitySquared(self):
"""
first eccentricity squared = (a**2-b**2)/a**2
"""
return (self.a ** 2 - self.b ** 2) / self.a ** 2
def getSecondEccentricity(self):
"""
second eccentricity = sqrt((a**2-b**2)/b**2)
"""
return math.sqrt((self.a ** 2 - self.b ** 2) / self.b ** 2)
def getSecondEccentricitySquared(self):
"""
second eccentricity squared (a**2-b**2)/b***2
"""
return (self.a ** 2 -self.b ** 2) / self.b ** 2
def getFlattering(self):
"""
flattering =(a-b)/a
"""
return (self.a - self.b) / self.a
def getInverseFlattering(self):
"""
inverse flattering =1/f = a/(a-b)
"""
return 1 / self.getFlattering()
def getNormal(self, lat):
"""
prime vertical, normal at latitude lat
"""
# return self.a**2/math.sqrt( \
# (self.a*math.cos(lat))**2 + \
# (self.b*math.sin(lat))**2)
return self.a / math.sqrt(\
1 - self.getFirstEccentricitySquared()*(math.sin(lat)) ** 2)
def getMeridional(self, lat):
"""
radius of curvature in the (north-south) meridian at latitude lat
"""
#-------------------------------- return (self.a*self.b)**2/math.pow(( \
#--------------------------------- (self.a*math.cos(lat))**2 + \
#--------------------------------- (self.b*math.sin(lat))**2), \
#-------------------------------------------------------- (3/2))
return (self.a * (1 - self.getFirstEccentricitySquared())) / (\
math.pow(1 - self.getFirstEccentricitySquared()*(math.sin(lat)) ** 2, 3 / 2))
def getPolarRadius(self):
"""
Distance from the Earth's center to the North (or South) Pole
"""
return self.b
def getEquatorialRadius(self):
"""
Distance from the Earth's center to the Equator
"""
return self.a
def getEllipsoidArea(self):
"""
Surface area S of an oblate ellipsoid
(generated by an ellipse rotating around its minor axis)
"""
# formula from:
# http://www.numericana.com/answer/geometry.htm#ovalarea
area=2 * math.pi * self.a ** 2 * (\
1 + (self.a / self.b) ** 2 * \
math.atan(self.getFirstEccentricity()) / \
(self.getFirstEccentricity()))
return area
def getEllipsoidVolume(self):
"""
Volume v of an oblate ellipsoid
(generated by an ellipse rotating around its minor axis)
"""
vol=(4 / 3) * math.pi * self.a ** 2 * self.b
return vol
def getAuthalicRadius(self):
"""
radius of a hypothetical perfect sphere which has
the same ellipsoid surface
"""
# formula from:
# http://en.wikipedia.org/wiki/Earth_radius#Notable_radii
return math.sqrt(self.getEllipsoidArea() / (math.pi * 4))
def getVolumetricRadius(self):
"""
radius of a hypothetical perfect sphere which has
the same ellipsoid volume
"""
# formula from:
# http://en.wikipedia.org/wiki/Earth_radius#Notable_radii
return math.pow(self.a**2*self.b, 1/3)
def getGeodeticRadius(self, lat):
"""
Distance from the Earth's center to a point on the
spheroid surface at geodetic latitude lat
"""
# formula from:
# http://en.wikipedia.org/wiki/Earth_radius#Notable_radii
gradius= math.sqrt(\
((self.a**2 * math.cos(lat))**2 + (self.b**2 * math.sin(lat))**2) / \
((self.a * math.cos(lat))** 2 + (self.b * math.sin(lat))** 2))
return gradius
def getParallelRadius(self, lat):
"""
radius of parallel at latitude lat
"""
return self.getNormal(lat) * math.cos(lat)
def getParallelArc(self, lat, fromLon, toLon):
"""
arc of parallel between two longitude (fromLong toLon) at laitutde lat
"""
return self.normal(lat) * math.cos(lat) * (toLon - fromLon)
def getMeridianArcDelambre(self,lat):
"""
arc of meridian from equator to lat
Delambre series
"""
e2, e4, e6, e8 = (math.pow(self.getFirstEccentricity(), i) for i in (2, 4, 6, 8))
k = self.a * (1 - e2)
A0 = k * (1 + 3 / 4 * e2 + 45 / 64 * e4 + \
175 / 256 * e6 + 11025 / 16384 * e8)
A2 = -(k / 2) * (3 / 4 * e2 + 15 / 16 * e4 + \
525 / 512 * e6 + 2205 / 2048 * e8)
A4 = k / 4 * (15 / 64 * e4 + 105 / 256 * e6 + 2205 / 4096 * e8)
A6 = -(k / 6) * (35 / 512 * e6 + 315 / 2048 * e8)
A8 = (k / 8) * (315 / 16348 * e8)
return A0*lat+A2*math.sin(2*lat)+A4*math.sin(4*lat)+ \
A6*math.sin(6*lat)+ A8*math.sin(8*lat)
def getMeridianArcBessel(self,lat):
"""
arc of meridian from equator to lat
Bessel, third flattering
"""
n=(self.a-self.b)/(self.a+self.b)
n2,n3,n4,n5=(math.pow(n, i) for i in range(2,6))
D0=1+(9/4*n**2)+(255/64*n4)
D2=(3/2*n)+(45/16*n3)+(525/128*n5)
D4=(15/16*n2)+(105/64*n4)
D6=(35/48*n3)+(315/256*n5)
return self.a*(1-n)**2*(1+n)*((D0*lat)-(D2*math.sin(2*lat))+\
(D4*math.sin(4*lat))-(D6*math.sin(6*lat)))
def getSferaLocale(self,lat):
"""raggio della sfera locale """
return math.sqrt(self.getMeridional(lat)*self.getNormal(lat))
if __name__ == "__main__":
'''qualche dato di prova'''
bessel=ellipsoid(6377397.155,f=0.003342773154,name='bessel')
wgs84 = ellipsoid(6378137.0, b=6356752.3142,name='wgs84')
grs80=ellipsoid(6378137.0,rf=298.257222101,name='grs80')
ell=[wgs84,bessel,grs80]
x=bessel.getFlattering()
lat=44
lat_rad=math.radians(lat)
print "Principali caratterisitche di alcuni ellissoidi"
print "Quando viene richiesto, la latitudine è uguale a {0} gradi".format(lat)
for e in ell:
print "========================================================="
print "Nome: {0}".format(e.name)
print "Flattering \'f\': {0}".format(e.getFlattering())
print "Inverse flatteirng \'rf\': {0}".format(e.getInverseFlattering())
print 'Eccentricità'.rjust(30,'-')
print "Eccentricità: {0}"\
.format(e.getFirstEccentricity())
print "Eccentricità al quadrato: {0}"\
.format(e.getFirstEccentricitySquared())
print "Eccentricità seconda: {0}"\
.format(e.getSecondEccentricity())
print "Eccentricità seconda al quadrato: {0}"\
.format(e.getSecondEccentricitySquared())
print 'Raggi'.rjust(30,'-')
print "Raggio equatoriale \'a\': {0}".format(e.a)
print "Raggio polare \'b\': {0}".format(e.b)
print "Raggio geodetico: {0:.15}"\
.format(e.getGeodeticRadius(lat_rad))
print "Gran normale alla latitudine di {0} gradi: {1}"\
.format(lat,e.getNormal(lat_rad))
print "Raggio di curvatura del meridiano alla latitudine di {0} gradi: {1}"\
.format(lat, e.getMeridional(lat_rad))
print "Raggio di curvatura del parallelo alla latitudine di {0} gradi: {1}"\
.format(lat, e.getParallelRadius(lat_rad))
print "Raggio di una sfera di area equivalente: {0}"\
.format(bessel.getAuthalicRadius())
print "Raggio di una sfera di volume equivalente: {0}"\
.format(bessel.getVolumetricRadius())
print "Raggio della sfera locale alla latitudine di {0} gradi: {1}"\
.format(lat,bessel.getSferaLocale(lat_rad))
print 'Archi di meridiano'.rjust(30,'-')
print "Arco di meridiano dall'equatore alla latitudine di {0} gradi (Bessel): {1}"\
.format(lat,e.getMeridianArcBessel(lat_rad))
print "Arco di meridiano dall'equatore alla latitudine di {0} gradi (Delambre): {1}"\
.format(lat,e.getMeridianArcDelambre(lat_rad))
0 commenti:
Posta un commento