Geodesy(Circular): Geodesy (circular) algorithms. Caution: This implement is largely untested. Tested on: Numworks 19.4.0 Numworks 20.3.2
"""Geodesyc: Geodesy (Circular) API. Note: Circular/Spherical calculations here. Ellipsoidal calculations can offer better accuracy (0.5% for Earth) (cf. Geodesy). Bearings are true N (not magnetic), clockwise. LatLongs are decimal (-ve S or W). Hint: on Earth, inputs of 4 decimals = 10 m nominal 5 decimals = 1 m Disclaimer: This implement is mainly untested. Refs https://www.movable-type.co.uk/scripts/latlong.html """ __version__ = "0.00" crid="""Geodesy V""" +__version__ +""" © 2023 SteveG1CMZ""" movable="""Some algorithms are adapted from https://www.movable-type.co.uk/scripts/latlong.html""" from math import * showme = True todeg=180/pi torad = pi/180 earthr=6371.000 #circular radius km kmtoNM=1.852 #nautical miles (exact) #algorithms def bearing1(lat1d,lon1d,lat2d,lon2d): """Initial bearing. Algorithm: from Movable. """ lat1= lat1d*torad lat2 = lat2d*torad lon1 = lon1d*torad lon2 = lon2d*torad deltalon=abs(lon2-lon1) #abs? y=sin(deltalon)*cos(lat2) x= cos(lat1)*sin(lat2) - sin(lat1)*cos(lat2)*cos(deltalon) theta = atan2(y,x) bearing= (theta*todeg + 360) % 360 return bearing #degrees def bearingback(bear1): return (bear1+180) % 360 def greatcircle(lat1d,lon1d,lat2d,lon2d): """great circle is Vincenty for great circles. Multiply by radius. Wiki says it's better conditioned than cos. Algorithm: from Wikipedia. """ lat1=lat1d*torad lat2=lat2d*torad lon1=lon1d*torad lon2=lon2d*torad coslat2=cos(lat2) deltalon=abs(lon2-lon1) cosdeltalon=cos(deltalon) T1=(coslat2*sin(deltalon))**2 T2=(cos(lat1)*sin(lat2) - sin(lat1)*cos(lat2)*cosdeltalon)**2 sq=T1+T2 den=sin(lat1)*sin(lat2) + cos(lat1)*coslat2*cosdeltalon delta = atan2(sqrt(sq),den) return delta #angular distance(radians, 0-pi) def gckm(lat1d,lon1d,lat2d,lon2d): dist=greatcircle(lat1d,lon1d,lat2d,lon2d) return (dist * earthr) def pt2(lat1d,lon1d,bearingd,delta): """Get pt2 given pt1 bearing and distance. delta=(angular) d/R, eg 100 km/earthr km Algorithm: from Movable. """ lat1 = lat1d*torad lon1 = lon1d*torad theta = bearing1*torad sinlat1=sin(lat1) #optimisations coslat1=cos(lat1) sindelta=sin(delta) cosdelta=cos(delta) lat2 = asin(sinlat1*cosdelta+ coslat1*sindelta*cos(theta)) lon2 = lon1 + atan2 (sin(theta)*sindelta*coslat1, cosdelta-sinlat1*sin(lat2)) lat2d=lat2*todeg lon2d = lon2*todeg lon2d = (lon2 + 540) % 360 - 180 #-180..180 return (lat2d, lon2d) def tosexa(degrees): mins = degrees-int(degrees) mins=60*mins secs= (mins-int(mins)) secs=60*secs return (int(degrees),int(mins), secs) def examples(): # a few examples #London 51.5072° N, 0.1276° W #Paris 48.8566 N, 2.3522 E km= gckm(51.5072,-0.1276,48.8566,2.3522) m= 1000* gckm(51.5072,-0.1276,48.8566,2.3522) NM=km/kmtoNM print(km,"km") print(m,"m") print(NM,"NM") print(gckm(40,0,60,0)) print(gckm(0,0,0,40)) print(gckm(0,90,0,-90)) km = gckm(51.5072,-0.1276,51.5072,-0.1276) print(km) bear1=bearing1(0,0,0,40) print(bear1,bearingback(bear1),"degrees") bear1=bearing1(0,0,90,0) print(bear1,bearingback(bear1),"degrees") if showme or __name__ == "__main__": print(crid) print(movable) examples()