Matplotlib Circle patch does not have smooth edges

user_9348121 0
Viewed 710

I'm trying to display Matplotlib patches using the Circle function on a map plot using cartopy geographical projections. Apparently this is supposed to give a smooth, near scale-free circular patch, however the edges are very polygonal. Strangely, CirclePolygon, the polygonal approximation counterpart of Circle, produces a smoother circle, albeit still not as smooth as I would like.

This is pretty much all the code as it pertains to adding the plot and the patches:

fig = plt.figure(figsize=(8,6))
img_extent = [340, 348, -35.5, -31]
ax = fig.add_subplot(1, 1, 1, projection = ccrs.Mollweide(), extent = img_extent)
patch_coords = [[342.5833, -34.5639],[343.4042, -34.3353],[343.8500, -33.8728],
                 [344.4917, -33.7636],[344.9250, -33.3108],[345.1333, -32.6811],
                   [344.9233, -32.1583]]
for pair in patch_coords:
    ax.add_patch(mpatches.Circle(xy = pair, radius = 0.5, 
                   color = 'r', alpha = 0.3, rasterized = None, 
                      transform = ccrs.Geodetic()))
ax.scatter(ra1, dec1, transform = ccrs.Geodetic(), rasterized = True, s = 1, 
             marker = ".", c = 'g', label = 'z < 0.025')
ax.scatter(ra2, dec2, transform = ccrs.Geodetic(), rasterized = True, s = 2, 
             marker = ".", c = 'b', label = '0.25 < z < 0.034')
ax.scatter(ra3, dec3, transform = ccrs.Geodetic(), rasterized = True, s = 0.75, 
             marker = ".", c = 'grey', label = '0.034 < z < 0.05')

Which produces this

Plot

I've tried looking through the available arguments but none seem to fix it. Is there a reason why it comes out like this and is there any way to make it smoother?

cartopy matplotlib python
2 Answers

I believe plotting Tissot's Indicatrices is more appropriate in your case. An Indicatrix represents a ground circle on a map projection. In many cases, the Indicatrices are rendered as ellipses as map projections do not always preserve shapes. The following is the working code that plots all the ground circles of radius = 55 km on the map projection that you desire. Read the comments in the code for some useful information.

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

fig = plt.figure(figsize=(12,8))

img_extent = [340, 348, -35.5, -31]
ax = fig.add_subplot(1, 1, 1, projection = ccrs.Mollweide(), extent = img_extent)

patch_coords = [[342.5833, -34.5639],[343.4042, -34.3353],[343.8500, -33.8728],
                 [344.4917, -33.7636],[344.9250, -33.3108],[345.1333, -32.6811],
                   [344.9233, -32.1583]]

for ix,pair in enumerate(patch_coords):
    # plot tissot indicatrix at each location 
    # n_samples = number of points forming indicatrix' perimeter
    # rad_km = 55 km. is about the angular distance 0.5 degree
    ax.tissot(rad_km=55, lons=np.array(patch_coords)[:,0][ix], \
              lats=np.array(patch_coords)[:,1][ix], n_samples=36, \
             facecolor='red', edgecolor='black', linewidth=0.15, alpha = 0.3)

gl = ax.gridlines(draw_labels=False, linewidth=1, color='blue', alpha=0.3, linestyle='--')
plt.show()

The resulting plot:

enter image description here

Edit

Since the first version of the code is not optimal. Code update is offered as follows:

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
fig = plt.figure(figsize=(12,8))

img_extent = [340, 348, -35.5, -31]
ax = fig.add_subplot(1, 1, 1, projection = ccrs.Mollweide(), extent = img_extent)

patch_coords = [[342.5833, -34.5639],[343.4042, -34.3353],[343.8500, -33.8728],
                 [344.4917, -33.7636],[344.9250, -33.3108],[345.1333, -32.6811],
                   [344.9233, -32.1583]]

for pair in patch_coords:
    # plot tissot indicatrix at each location 
    # n_samples = number of points forming indicatrix' perimeter
    # rad_km = 55 km. is about the angular distance 0.5 degree at equator
    ax.tissot(rad_km=55, lons=pair[0], lats=pair[1], n_samples=36, \
             facecolor='red', edgecolor='black', linewidth=0.15, alpha = 0.3)

gl = ax.gridlines(draw_labels=False, linewidth=1, color='blue', alpha=0.3, linestyle='--')
plt.show()
user_21774130

I believe that Cartopy does line projections with an arbitrary fixed accuracy, rather than a dynamic line-split calculation.

See e.g. :

I also think work is ongoing right now to address that.

In the meantime, to solve specific problems you can hack the CRS.threshold property, as explained here : https://github.com/SciTools/cartopy/issues/8
That is, you can make it use finer steps by reprogramming the fixed value.
I think this would also fix this circle-drawing problem, though I'm not 100%

user_26150500
Related