TM5
TM5

Reputation: 302

How to plot a filled polygon on a map in cartopy

Edited, adding suggestion from an answer

I have a list of vertices in lat/lon that define corners of a polygon on a map. I would like to draw that polygon on a map using cartopy, where the edges are great circles. I've tried following the examples at https://scitools.org.uk/cartopy/docs/v0.5/matplotlib/introductory_examples/02.polygon.html, but I can't get it to work. Here's what I have tried so far:

import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import matplotlib.patches as mpatches

map_proj = ccrs.Orthographic(central_latitude=0.0, central_longitude=80.0)
ax = plt.axes(projection=map_proj)

ax.set_global() # added following an answer to my question
ax.gridlines()

ax.coastlines(linewidth=0.5, color='k', resolution='50m')

lat_corners = np.array([-20.,  0., 50., 30.])
lon_corners = np.array([ 20., 90., 90., 30.]) + 15.0 # offset from gridline for clarity

poly_corners = np.zeros((len(lat_corners), 2), np.float64)
poly_corners[:,0] = lon_corners
poly_corners[:,1] = lat_corners

poly = mpatches.Polygon(poly_corners, closed=True, ec='r', fill=False, lw=1, fc=None, transform=ccrs.Geodetic())
ax.add_patch(poly)

Output image

Notice that the lines are not great circles, and there seem to be more than four vertices. I feel like this is such a simple thing to do there must be a way, but I can't figure that out from the cartopy documentation.

Upvotes: 3

Views: 4325

Answers (3)

swatchai
swatchai

Reputation: 18772

Here is a runnable code and output. Credits go to: the asker, ImportanceOfBeingErnest, and ajdawson.

import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import matplotlib.patches as mpatches

map_proj = ccrs.Orthographic(central_latitude=0.0, central_longitude=80.0)

map_proj._threshold /= 100.  # the default values is bad, users need to set them manually

ax = plt.axes(projection=map_proj)

ax.set_global() # added following an answer to my question
ax.gridlines()

ax.coastlines(linewidth=0.5, color='k', resolution='50m')

lat_corners = np.array([-20.,  0., 50., 30.])
lon_corners = np.array([ 20., 90., 90., 30.]) + 15.0 # offset from gridline for clarity

poly_corners = np.zeros((len(lat_corners), 2), np.float64)
poly_corners[:,0] = lon_corners
poly_corners[:,1] = lat_corners

poly = mpatches.Polygon(poly_corners, closed=True, ec='r', fill=True, lw=1, fc="yellow", transform=ccrs.Geodetic())
ax.add_patch(poly)

plt.show()

Output:

enter image description here

Upvotes: 2

ajdawson
ajdawson

Reputation: 3333

I think this is probably because Cartopy's default transform resolution is too low for this projection. You can work around this by forcing a higher resolution:

map_proj = ccrs.Orthographic(central_latitude=0.0, central_longitude=80.0)
map_proj._threshold /= 100.
...

This gives nice curved great circle arcs.

Upvotes: 3

ImportanceOfBeingErnest
ImportanceOfBeingErnest

Reputation: 339112

Mind, that the example uses

ax.set_global()

Upvotes: 1

Related Questions