Skip to content

Instantly share code, notes, and snippets.

@tianyeeee
Forked from arthur-e/geodesic.py
Created April 12, 2021 09:43
Show Gist options
  • Save tianyeeee/9f2109490fb84b71fcdd2d62bb56174e to your computer and use it in GitHub Desktop.
Save tianyeeee/9f2109490fb84b71fcdd2d62bb56174e to your computer and use it in GitHub Desktop.
A basic example of creating a geotrajectory visualization in KML, with a great circle drawn between two points of interest. Requires geodesic.py, a script that provides a function "tweensegs."
import math
def tweensegs(longitude1, latitude1, longitude2, latitude2, num_of_segments):
ptlon1 = longitude1
ptlat1 = latitude1
ptlon2 = longitude2
ptlat2 = latitude2
numberofsegments = num_of_segments
onelessthansegments = numberofsegments - 1
fractionalincrement = (1.0/onelessthansegments)
ptlon1_radians = math.radians(ptlon1)
ptlat1_radians = math.radians(ptlat1)
ptlon2_radians = math.radians(ptlon2)
ptlat2_radians = math.radians(ptlat2)
distance_radians=2*math.asin(math.sqrt(math.pow((math.sin((ptlat1_radians-ptlat2_radians)/2)), 2) + math.cos(ptlat1_radians)*math.cos(ptlat2_radians)*math.pow((math.sin((ptlon1_radians-ptlon2_radians)/2)), 2)))
# 6371.009 represents the mean radius of the earth
# shortest path distance
distance_km = 6371.009 * distance_radians
mylats = []
mylons = []
# Write the starting coordinates
mylats.append([])
mylons.append([])
mylats[0] = ptlat1
mylons[0] = ptlon1
f = fractionalincrement
icounter = 1
while (icounter < onelessthansegments):
icountmin1 = icounter - 1
mylats.append([])
mylons.append([])
# f is expressed as a fraction along the route from point 1 to point 2
A=math.sin((1-f)*distance_radians)/math.sin(distance_radians)
B=math.sin(f*distance_radians)/math.sin(distance_radians)
x = A*math.cos(ptlat1_radians)*math.cos(ptlon1_radians) + B*math.cos(ptlat2_radians)*math.cos(ptlon2_radians)
y = A*math.cos(ptlat1_radians)*math.sin(ptlon1_radians) + B*math.cos(ptlat2_radians)*math.sin(ptlon2_radians)
z = A*math.sin(ptlat1_radians) + B*math.sin(ptlat2_radians)
newlat=math.atan2(z, math.sqrt(math.pow(x, 2)+math.pow(y, 2)))
newlon=math.atan2(y, x)
newlat_degrees = math.degrees(newlat)
newlon_degrees = math.degrees(newlon)
mylats[icounter] = newlat_degrees
mylons[icounter] = newlon_degrees
icounter += 1
f = f + fractionalincrement
# Write the ending coordinates
mylats.append([])
mylons.append([])
mylats[onelessthansegments] = ptlat2
mylons[onelessthansegments] = ptlon2
# Now, the array mylats[] and mylons[] have the coordinate pairs for intermediate points along the geodesic
# My mylat[0], mylat[0] and mylat[num_of_segments-1], mylat[num_of_segments-1] are the geodesic end points
buff = ''
# Write a kml of the results
zipcounter = 0
kmlheader = "<?xml version=\"1.0\" encoding=\"UTF-8\"?><kml xmlns=\"http://www.opengis.net/kml/2.2\"><Document><name>LineString.kml</name><open>1</open><Placemark><name>unextruded</name><LineString><extrude>1</extrude><tessellate>1</tessellate><coordinates>"
buff += kmlheader
while (zipcounter < numberofsegments):
outputstuff = repr(mylons[zipcounter]) + "," + repr(mylats[zipcounter]) + ",0 "
buff += outputstuff
zipcounter += 1
kmlfooter = "</coordinates></LineString></Placemark></Document></kml>"
buff += kmlfooter
return buff
import ogr, pykml
from pykml.parser import fromstring as from_kml_string
from pykml.factory import KML_ElementMaker as KML
from lxml import etree
from geodesic import tweensegs
# For geodesic.py and the function "tweensegs" see this article:
# http://gis.stackexchange.com/questions/47/what-tools-in-python-are-available-for-doing-great-circle-distance-line-creati/316#316
# These coordinates must be strings to keep precision
PARIS = ['2.35', '48.86', 'Paris, France']
ROME = ['12.50', '41.90', 'Rome, Italy']
NYC = ['-73.94', '40.66', 'New York City, U.S.A.']
JERUSALEM = ['35.22', '31.78', 'Jerusalem']
PT_A = PARIS
PT_B = JERUSALEM
pt_a = ogr.CreateGeometryFromWkt('POINT(%s %s)' % tuple(PT_A[0:2]))
pt_b = ogr.CreateGeometryFromWkt('POINT(%s %s)' % tuple(PT_B[0:2]))
# Buffer the points, an ellipse with radius 0.1 degrees, with 20 points per arc
buff_a = pt_a.Buffer(10.0, 20)
buff_b = pt_b.Buffer(15.0, 20)
features = [
buff_a,
buff_b
]
placemarks = [
KML.Placemark(
KML.name(PT_A[2]),
KML.Point(
KML.coordinates('%s, %s, 0' % tuple(PT_A[0:2]))
),
KML.Style(
KML.IconStyle(
KML.Icon(
KML.href('http://maps.google.com/mapfiles/kml/paddle/A.png')
),
KML.scale(4)
)
)
),
KML.Placemark(
KML.name(PT_B[2]),
KML.Point(
KML.coordinates('%s, %s, 0' % tuple(PT_B[0:2]))
),
KML.Style(
KML.IconStyle(
KML.Icon(
KML.href('http://maps.google.com/mapfiles/kml/paddle/B.png')
),
KML.scale(4)
)
)
),
KML.Placemark(
KML.Style(
KML.LineStyle(
KML.color('ffffffff'),
KML.width(3)
)
),
KML.LineString(
KML.coordinates(str(from_kml_string(tweensegs(float(PT_A[0]), float(PT_A[1]), float(PT_B[0]), float(PT_B[1]),
20)).xpath('kml:Document/kml:Placemark/kml:LineString/kml:coordinates',
namespaces={'kml': "http://www.opengis.net/kml/2.2"})[0]).strip())
)
)
]
for feature in features:
# Parse each feature as a KML string and append the object to a sequence
placemarks.append(KML.Placemark(
KML.Style(
KML.LineStyle(
KML.color('ff0000ff'),
KML.width(3)
),
KML.PolyStyle(
KML.color('770000ff')
)
),
from_kml_string(feature.ExportToKML())
))
document = KML.Document(*placemarks)
with open('GeoTrajectory.kml', 'wb') as kml_doc:
kml_doc.write(etree.tostring(document, pretty_print=True))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment