¿Qué herramientas en Python están disponibles para hacer una gran distancia de círculo + creación de línea?

20

Necesito usar Python para crear una gran distancia circular, tanto un número como, preferiblemente, algún tipo de 'curva' que pueda usar para dibujar en un mapa del lado del cliente. No me importa el formato de la curva, ya sea WKT o un conjunto de pares de coordenadas, pero solo quiero obtener los datos.

¿Qué herramientas hay por ahí? ¿Qué debo usar?

Christopher Schmidt
fuente

Respuestas:

8

Las respuestas proporcionadas por otros son un poco más elegantes, pero aquí hay un poco de Python ultra simple, un tanto no pitónico que proporciona lo básico. La función toma dos pares de coordenadas y un número de segmentos especificado por el usuario. Produce un conjunto de puntos intermedios a lo largo de un gran camino circular. Salida: texto listo para escribir como KML. Advertencias: El código no considera antípodas, y asume una tierra esférica.

Código de Alan Glennon http://enj.com Julio de 2010 (el autor coloca este código en el dominio público. Utilícelo bajo su propio riesgo).

-

def tweensegs (longitud1, latitud1, longitud2, latitud2, número_de_segmentos):

import math

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

# 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>"
print kmlheader
while (zipcounter < numberofsegments):
        outputstuff = repr(mylons[zipcounter]) + "," + repr(mylats[zipcounter]) + ",0 "
        print outputstuff
        zipcounter += 1
kmlfooter = "</coordinates></LineString></Placemark></Document></kml>"
print kmlfooter
Glennon
fuente
8

GeographicLib tiene una interfaz de python :

Esto puede computar geodésicas en un elipsoide (establecer el aplanamiento a cero para obtener grandes círculos) y puede generar puntos intermedios en una geodésica (ver los comandos "Línea" en la muestra).

Aquí le mostramos cómo imprimir puntos en la línea geodésica de JFK al aeropuerto de Changi (Singapur):

from geographiclib.geodesic import Geodesic
geod = Geodesic.WGS84

g = geod.Inverse(40.6, -73.8, 1.4, 104)
l = geod.Line(g['lat1'], g['lon1'], g['azi1'])
num = 15  # 15 intermediate steps

for i in range(num+1):
    pos = l.Position(i * g['s12'] / num)
    print(pos['lat2'], pos['lon2'])

->
(40.60, -73.8)
(49.78, -72.99)
(58.95, -71.81)
(68.09, -69.76)
(77.15, -65.01)
(85.76, -40.31)
(83.77, 80.76)
(74.92, 94.85)
...
cffk
fuente
El puerto de Python de GeographicLib ahora está disponible en pypi.python.org/pypi/geographiclib
cffk
Vea también este documento: CFF Karney, Algorithms for Geodesics, J. Geod, DOI: dx.doi.org/10.1007/s00190-012-0578-z
cffk
7

pyproj tiene la función Geod.npts que devolverá una serie de puntos a lo largo del camino. Tenga en cuenta que no incluye los puntos terminales en la matriz, por lo que debe tenerlos en cuenta:

import pyproj
# calculate distance between points
g = pyproj.Geod(ellps='WGS84')
(az12, az21, dist) = g.inv(startlong, startlat, endlong, endlat)

# calculate line string along path with segments <= 1 km
lonlats = g.npts(startlong, startlat, endlong, endlat,
                 1 + int(dist / 1000))

# npts doesn't include start/end points, so prepend/append them
lonlats.insert(0, (startlong, startlat))
lonlats.append((endlong, endlat))
scruss
fuente
¡Gracias! Solución provista por una biblioteca conocida y utilizada masivamente aquí :)
tdihp
-2

No he usado este paquete, pero parece interesante y una posible solución: http://trac.gispython.org/lab/wiki/Shapely

Hugo Estrada
fuente
77
AFAIK, bien proporcionado no hace cálculos esferoidales:Shapely is a Python package for set-theoretic analysis and manipulation of **planar** features
fmark