# Getting geometries of hurricane's 'cone of uncertainty' using shapely?

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 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 = [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))])
``````

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!

• try looking for convex-hull or alpha-shapes – Ian Turton Jun 24 '19 at 7:30
• What you need is a variable width buffer that connects the tangents of subsequent circles, see twitter.com/geospacedman/status/978172211579228162 for how QGIS does it. This is not built into Shapely as far as I know, would be super useful though if you fancy a PR. – bugmenot123 Jun 24 '19 at 8:19

Shapely geometries have a `convex_hull` method.

Should be as simple as `polygons.convex_hull`, but it will work with any Shapely geometry.

A note on cyclones as a domain: you should use the input cyclone positions as input rather than an interpolated curve: weather forecasts are typically made for moments in time, often 3, 6 or 12 hours apart, and the position in-between is uncertain simply because it is left uncalculated. The convex hull (a special kind of alpha shape) will encompass the spaces in-between the foretasted locations, exactly like in your sample images.

Also be careful with the antimeridian...

Edit: on second thought, you probably want a concave hull, or else to generate convex hulls sequentially, starting with the first pair of error shapes, then with the i+1 and i+2 pair, until complete. Then you union this set of pair-wise convex hulls together. If you do a simple convex hull, then your overall shape will be, well, convex rather than somewhat concave. But a naive concave hull may well be too "tight" and cause intrusions into the path that you don't want.

To illustrate (pseudo-code):

``````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)
``````
• Convex Hull would ruin the concave parts at the top of OP's image though. – bugmenot123 Jun 24 '19 at 8:10
• Yes, I just realised. It's still possible to use a convex hull, but it needs to be a progressive operation, like a list `reduce`. See edit. – alphabetasoup Jun 24 '19 at 8:14

If you need a polygon like in the image below, replace the last lines of your code by the following:

``````#### 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
``````

Convex hulls: Final result: There is an implementation of `VariableWidthBuffer` in the JTS Lab here. It uses a union of circles around each line vertex and "prisms" around each line segment. That could be a basis for a Python implementation.

This will make it into JTS sometime soon. Then perhaps into GEOS, where it can be exposed by Shapely.

• Actually reviewing the Lab code I realized that the approach there is too simplistic and does not produce quality buffer outlines. I have a new approach coded up and will be committing that to JTS soon . (Warning: much more math involved!) – dr_jts Nov 19 '19 at 20:22
• An improved `VariableBuffer` implementation is now available in JTS here – dr_jts Nov 20 '19 at 20:30
• It shouldn't be too hard to port the JTS VariableBuffer code to Python, using Shapely for the geometric heavy lifting. – dr_jts Nov 20 '19 at 20:57