Determinando el punto de un “Gran Círculo” con la menor distancia a un punto exterior

Cuando se trabaja con coordenadas proyectadas referidas a pequeñas distancias, la separación más corta entre dos puntos (expresada mediante una línea adecuadamente densificada) probablemente no difiera visualmente de una recta. Sin embargo, cuando éstas se incrementan, como la representación de la superficie de la tierra corresponde a la de un elipsoide, la menor distancia entre dos puntos, en la mayoría de las proyecciones que comunmente usamos, debería visualizarse como una curva; no como una recta. Estas curvas corresponden a lo que denominamos “Grandes Círculos”.

Los “Grandes Círculos” se visualizan como rectas en las proyecciones azimutales equidistantes o gnomicas por lo que estas se emplean para crearlos.

En la imagen siguiente, se observa (en naranja) un “Gran Círculo” creado entre el punto (7,51) y uno arbitrario sobre los Estados Unidos; usado el primero como base para crear una proyección azimutal equidistante (aeqd) personalizada:

+proj=aeqd +lat_0=51 +lon_0=7 +x_0=0 +y_0=0 +a=6371000 +b=6371000 +units=m +no_defs

La geometría tipo línea se densificó en QGIS con 100 puntos intermedios antes de desplegarla con un EPSG:4326 (long/lat WGS 84).

great_circle1

La línea en verde representa una ruta arbitraria trazada sobre el mapa y se desea determinar los puntos en el “Gran Círculo” que representan la distancia más corta hasta los puntos marcados sobre la referida línea verde. El código a continuación realiza el proceso imprimiendo en la Python Console las coordenadas de tales puntos y los segmentos de recta que los unen (que no tienen nada que ver con sus respectivos “Grandes Círculos”).

mapcanvas = iface.mapCanvas()

layers = mapcanvas.layers()

#for point layer
feats_points = [  feat for feat in layers[0].getFeatures()  ] 

#for great circle
feat_gc = layers[1].getFeatures().next()

lines  = []
red_points = []

for feat in feats_points:
    lines.append([feat.geometry().asPoint(), feat_gc.geometry().closestSegmentWithContext(feat.geometry().asPoint())[1]])
    red_points.append(feat_gc.geometry().closestSegmentWithContext(feat.geometry().asPoint())[1])

print QgsGeometry.fromMultiPoint(red_points).exportToWkt()
print QgsGeometry.fromMultiPolyline(lines).exportToWkt()

Después de ejecutado el código anterior en la Python Console de QGIS se tiene:

great_circle2

Con la ayuda del plugin QuickWKT los WKT resultantes se visualizaron de manera expedita en el Map Canvas de QGIS y representan, respectivamente, los puntos en rojo y los segmentos de recta que se observan en la imagen superior.

Esta entrada fue publicada en PyQGIS, SIG, Software Libre. Guarda el enlace permanente.

Responder

Por favor, inicia sesión con uno de estos métodos para publicar tu comentario:

Logo de WordPress.com

Estás comentando usando tu cuenta de WordPress.com. Cerrar sesión / Cambiar )

Imagen de Twitter

Estás comentando usando tu cuenta de Twitter. Cerrar sesión / Cambiar )

Foto de Facebook

Estás comentando usando tu cuenta de Facebook. Cerrar sesión / Cambiar )

Google+ photo

Estás comentando usando tu cuenta de Google+. Cerrar sesión / Cambiar )

Conectando a %s