8

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!

2
  • try looking for convex-hull or alpha-shapes
    – Ian Turton
    Commented Jun 24, 2019 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. Commented Jun 24, 2019 at 8:19

3 Answers 3

11

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)
2
  • Convex Hull would ruin the concave parts at the top of OP's image though. Commented Jun 24, 2019 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. Commented Jun 24, 2019 at 8:14
7

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

2

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.

3
  • 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
    Commented Nov 19, 2019 at 20:22
  • 1
    An improved VariableBuffer implementation is now available in JTS here
    – dr_jts
    Commented Nov 20, 2019 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
    Commented Nov 20, 2019 at 20:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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