7

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:

1 2

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:

pic

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 at 7:30
  • 1
    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 at 8:19
7

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 at 8:10
  • 1
    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 at 8:14
5

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
polygons = cascaded_union(convex_hulls)

Convex hulls:

enter image description here

Final result:

enter image description here

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.