AxelEckenberger
AxelEckenberger

Reputation: 16926

Plot a circle on a cartopy mercator projection

For a project I need to create a visualization that draws a circle around some locations on a map. The visualization used Cartopy v.0.18.0 to render the map. It uses the GoogleTiles class to fetch and display the tiles in the relevant region, and the add_patch(Patch.Circle(..., transform=ccrs.PlateCarree())) method to draw the circle.

tiles = GoogleTiles()
fig = plt.figure(figsize=(15,15))
ax = fig.add_subplot(1, 1, 1, projection=tiles.crs)

ax.set_extent((-121.8,-122.55,37.25,37.85))

ax.add_image(tiles, 11)

ax.add_patch(Patch.Circle(xy=[-122.4015173428571, 37.78774634285715], radius = 0.021709041989311614 + 0.005, alpha=0.3, zorder=30, transform=ccrs.PlateCarree()))

plt.show()

However, although I tried several transform objects I either got a ellipse instead of a circle (e.g. using ccrs.PlateCarree()) or no circle at all (e.g. using ccrs.Mercator()).

I found several different solutions online (e.g. Drawing Circles with cartopy in orthographic projection), however, these were not for the Mercator projection and I sadly lack the projection/transformation knowledge to adapt these to my problem.

Circle on Mercator projection

The only way I was able to produce a circular patch, was when I set the projection parameter on fig.add_subplot to ccrs.PlateCarree(). This, however, distorts the map and the labels become blured, so this is sadly not an acceptable solution.

As the project is due soon, a speedy reply would be much appreciated.

Upvotes: 2

Views: 1488

Answers (1)

AxelEckenberger
AxelEckenberger

Reputation: 16926

Thanks @swatchai this was the missing hint, so for those intested the code looks like this right now, and it does work! Hooray!

tiles = GoogleTiles()
fig = plt.figure(figsize=(15,15))
ax = fig.add_subplot(1, 1, 1, projection=tiles.crs)

ax.set_extent((-121.8,-122.55,37.25,37.85))

ax.add_image(tiles, 11)

# The diameter is in degrees in EPSG:4326 coordinates therefore, the degrees have 
# to be converted to km. At 37N the degree latitude is 11.0977 km.
ax.tissot(rad_km=(0.021709041989311614 + 0.005) * 11.0977, lons=[-122.4015], lats=[37.7877], alpha=0.3)

plt.show()

When executing the above code the following warning is thrown but it has visible effect on the result:

/opt/conda/lib/python3.8/site-packages/cartopy/mpl/geoaxes.py:761: UserWarning: Approximating coordinate system <cartopy._crs.Geodetic object at 0x7fa4c7529770> with the PlateCarree projection.
  warnings.warn('Approximating coordinate system {!r} with the '

Circle on Mercator map

So thanks again @swatchai you saved my day!

Upvotes: 2

Related Questions