[FIXED] How to draw the radius of a circle within a Cartopy projection

Issue

I am trying to draw the radius of a circle on a Cartopy projection through one point. And I couldn't find anything related to this.

Here is what I have so far:

fig = plt.figure(figsize=(10, 6))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.Mercator())

ax.set_extent([18, 28, 59.5, 64.1], crs=ccrs.PlateCarree())
ax.coastlines(linewidth=.5)

# Add the radar distance circle
lon_ika = 23.076
lat_ika = 61.76
radius = 250
n_samples = 80

circles = Polygon(Geodesic().circle(lon_ika, lat_ika, radius*1000., n_samples=n_samples))
feature = cfeature.ShapelyFeature(circles, ccrs.PlateCarree(), fc='None', ec="black", lw=1, linestyle="-")
linestyle="--")
circle = ax.add_feature(feature)

# Adding red dot and name of the radar station to the plot
ax.plot(lon_ika, lat_ika, "o", c='r', transform=ccrs.PlateCarree(), markersize=6, label="Ikaalinen")

# Adding red cross and name of IOP location to the plot
lon_hyy = 24.3
lat_hyy = 61.83
ax.plot(lon_hyy, lat_hyy, "x", c='r', transform=ccrs.PlateCarree(), markersize=6, label="Hyytiälä")

# Add labels
plt.legend(loc='upper left', fontsize=12, framealpha=1, edgecolor='black')

plt.show()

I haven't found anything related on the web so far :/

Solution

This code should give the required radius line that passes the X mark.

The code:

from shapely.geometry import Polygon
from cartopy.geodesic import Geodesic  # from geographiclib
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy
import cartopy.feature as cfeature

fig = plt.figure(figsize=(10, 6))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.Mercator())

# Avoid error AttributeError: 'GeoAxesSubplot' object has no attribute '_autoscaleXon'
ax._autoscaleXon = False
ax._autoscaleYon = False

ax.set_extent([18, 28, 59.5, 64.1], crs=ccrs.PlateCarree())
ax.coastlines(linewidth=.5)

# Add the radar distance circle
lon_ika = 23.076
lat_ika = 61.76
ika = [lon_ika, lat_ika]
radius = 250     # km
n_samples = 80

circles = Polygon(Geodesic().circle(lon_ika, lat_ika, radius*1000., n_samples=n_samples))
feature = cfeature.ShapelyFeature([circles], ccrs.PlateCarree(), fc='None', ec="black", lw=1, linestyle="-")

circle = ax.add_feature(feature)

# Adding red dot and name of the radar station to the plot
ax.plot(lon_ika, lat_ika, "o", c='r', transform=ccrs.PlateCarree(), markersize=6, label="Ikaalinen", zorder=30)

# Adding red cross and name of IOP location to the plot
lon_hyy = 24.3
lat_hyy = 61.83
hyy = [lon_hyy, lat_hyy]

# Get (geod_distance, forward and backward azimuth) between 2 points
dist_m, fw_azim, bw_azim  = Geodesic().inverse(ika, hyy).T
# Get (long, lat, forward_azimuth) of target point using direct problem solver
px_lon, px_lat, fwx_azim = Geodesic().direct(ika, fw_azim, radius*1000).T

ax.plot(lon_hyy, lat_hyy, "x", c='r', transform=ccrs.PlateCarree(), markersize=6, label="Hyytiälä", zorder=10)

# Plot the target point on the circle's perimeter
ax.plot(px_lon, px_lat, "x", c='g', transform=ccrs.PlateCarree(), markersize=12, label="Target")

# Plot great-circle arc from circle center to the target point
ax.plot([ika[0], px_lon], [ika[1], px_lat], '.-', color='blue',  transform=ccrs.Geodetic(), zorder=5 )

gl = ax.gridlines(draw_labels=True)

#ax.set_aspect(1)

# Add labels
plt.legend(loc='upper left', fontsize=12, framealpha=1, edgecolor='black')

plt.show()

The output plot:

Answered By - swatchai