Obtención de geometrías del "cono de incertidumbre" del huracán mediante shapely?

Question

Estoy tratando de trazar mapas de previsión de huracanes usando python. Tengo varias posiciones de previsión derivadas de los avisos oficiales, y las interpolo en una curva suavizada, luego dibujo un polígono de "cono de incertidumbre" basado en la curva. Ejemplo:

Básicamente, el "cono de incertidumbre" es la huella de un círculo que se mueve y se amplía. He probado muchos enfoques, pero ninguno de ellos es lo suficientemente bueno. Mi enfoque actual es producir ~100 círculos basados en la curva interpolada, y hacer un polígono compuesto usando cascaded_union método en shapely .

import numpy as np
from shapely.geometry import MultiPolygon
from shapely.ops import cascaded_union
from scipy.interpolate import interp1d
# x, y: coords of forecast position
y = [18.3, 19.2, 20.0, 20.4, 20.7, 21.3, 21.6, 21.5, 20.8, 20.8, 21.5]
x = [111.3, 111.2, 110.9, 110.5, 110.2, 110.5, 110.0, 109.2, 109.4, 110.3, 111.8]
# r: radius of uncertainty
r = [0, 0.05, 0.1, 0.15, 0.2, 0.3, 0.4, 0.5, 0.5, 0.5, 0.5]
hours = [0, 6, 12, 18, 24, 36, 48, 60, 72, 96, 120]
# interpolate
points_num = 100
interp_hours = np.linspace(min(hours), max(hours), points_num)
x = interp1d(hours, x, kind='cubic')(interp_hours)
y = interp1d(hours, y, kind='cubic')(interp_hours)
r = interp1d(hours, r, kind='linear')(interp_hours)
# make polygon
thetas = np.linspace(0, 2 * np.pi, 360)
polygon_x = x[:,None] + r[:,None] * np.sin(thetas)
polygon_y = y[:,None] + r[:,None] * np.cos(thetas)
polygons = MultiPolygon([Polygon(i) for i in np.dstack((polygon_x, polygon_y))])
polygons = cascaded_union(polygons).buffer(0)

Pero se ve mal cerca del punto de partida:

Aumentar el número de círculos sólo puede resolver parcialmente el problema y lleva más tiempo. Así que me pregunto si hay una hermosa, eficiente y pitonisa ¿se puede hacer el "cono de incertidumbre"? Hay que tener en cuenta que los huracanes pueden cambiar sus direcciones bruscamente, ¡e incluso estar estacionados!

Answer 1

Las geometrías con forma tienen un convex_hull método .

Debería ser tan sencillo como polygons.convex_hull pero funcionará con cualquier geometría Shapely.

Una nota sobre los ciclones como dominio: hay que utilizar las posiciones de los ciclones de entrada en lugar de una curva interpolada: las previsiones meteorológicas se hacen normalmente para momentos en el tiempo, a menudo con 3, 6 o 12 horas de diferencia, y la posición intermedia es incierta simplemente porque se deja sin calcular. El casco convexo (un tipo especial de forma alfa ) englobará los espacios entre las ubicaciones de los predios, exactamente como en sus imágenes de muestra.

También hay que tener cuidado con el antimeridiano ...

Edición: pensándolo bien, probablemente quieras un casco cóncavo, o bien generar cascos convexos secuencialmente, empezando con el primer par de formas de error, luego con el par i+1 e i+2, hasta completarlo. Entonces se une este conjunto de cascos convexos por pares. Si haces un casco convexo simple, entonces tu forma general será convexo en lugar de algo cóncavo. Pero un casco cóncavo ingenuo puede ser demasiado "apretado" y causar intrusiones en la trayectoria que no se desea.

Para ilustrar (pseudocódigo):

shapes = [a, b, c, d] # Ordered list of shapely geometries
parts = []
for first, second in zip(shapes, shapes[1:]):
    parts.append(union(first, second).convex_hull)
union(*parts)

Answer 2

Si necesitas un polígono como el de la imagen de abajo, sustituye las últimas líneas de tu código por lo siguiente:

#### firstly, import Polygon class ####
from shapely.geometry import MultiPolygon, Polygon
.
.
.
# make polygon
thetas = np.linspace(0, 2 * np.pi, 360)
polygon_x = x[:,None] + r[:,None] * np.sin(thetas)
polygon_y = y[:,None] + r[:,None] * np.cos(thetas)

# circles
ps = [Polygon(i) for i in np.dstack((polygon_x, polygon_y))]

# list of convex hulls of subsequent circles
n = range(len(ps)-1)
convex_hulls = [MultiPolygon([ps[i], ps[i+1]]).convex_hull for i in n]

# Final polygon
polygons = cascaded_union(convex_hulls)

Cascos convexos:

Resultado final:

Answer 3

Existe una aplicación de VariableWidthBuffer en el laboratorio del STC aquí . Utiliza una unión de círculos alrededor de cada vértice de línea y "prismas" alrededor de cada segmento de línea. Eso podría ser una base para una implementación en Python.

Esto llegará a la STC en algún momento. Luego, tal vez en GEOS, donde puede ser expuesto por Shapely.