Search code examples
pythonmatplotlibvisualizationgeospatialcartopy

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.


Solution

  • 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!