I'm trying to plot hurricane forecast maps using python. I have several forecast positions derived from official advisories, and interpolate them into a smoothed curve, then draw a polygon of 'cone of uncertainty' based on the curve. Example:
Basically, 'cone of uncertainty' is the footprint of a moving and enlarging circle. I've tried many approaches but none of them is good enough. My current approach is to produce ~100 circles based on the interpolated curve, and make a compound polygon using cascaded_union
method in 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)
But it looks nasty near the starting point:
Increasing the number of circles can only partly address the problem and it takes more time. So I wonder if there is a beautiful, efficient and pythonic way to make the 'cone of uncertainty'? Note that hurricanes may change its directions abruptly, and even be stationary!