Drawing circles with Cartopy in NorthPolarStereo projection

Question

I would like to draw circles in Cartopy in NorthPolarStereo projections, providing the center and radius in lat,lon units.

Similar and excellent questions and answers are available for Basemap here, and for Cartopy in Ortographic projection here. However, I would like to use the NorthPolarStereo in Cartopy. Trying the latter approach, just changing the projection makes the circle fixed in the North Pole, ignoring the coordinates you give for its center.

Any ideas on how to draw circles in Cartopy using NorthPolarStereo projection, providing its center and radius in lat,lon untis?

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

# example: draw circle with 45 degree radius around the North pole
lat = 72 
lon = 100
r = 20

# Define the projection used to display the circle:
proj = ccrs.NorthPolarStereo(central_longitude=lon)


def compute_radius(ortho, radius_degrees):
    phi1 = lat + radius_degrees if lat <= 0 else lat - radius_degrees
    _, y1 = ortho.transform_point(lon, phi1, ccrs.PlateCarree())
    return abs(y1)

# Compute the required radius in projection native coordinates:
r_ortho = compute_radius(proj, r)

# We can now compute the correct plot extents to have padding in degrees:
pad_radius = compute_radius(proj, r + 5)

# define image properties
width = 800
height = 800
dpi = 96
resolution = '50m'

# create figure
fig = plt.figure(figsize=(width / dpi, height / dpi), dpi=dpi)
ax = fig.add_subplot(1, 1, 1, projection=proj)

ax.set_xlim([-pad_radius, pad_radius])
ax.set_ylim([-pad_radius, pad_radius])
ax.imshow(np.tile(np.array([[cfeature.COLORS['water'] * 255]], dtype=np.uint8), [2, 2, 1]), origin='upper', transform=ccrs.PlateCarree(), extent=[-180, 180, -180, 180])
ax.add_feature(cfeature.NaturalEarthFeature('physical', 'land', resolution, edgecolor='black', facecolor=cfeature.COLORS['land']))
ax.add_patch(mpatches.Circle(xy=[lon, lat], radius=r_ortho, color='red', alpha=0.3, transform=proj, zorder=30))

plt.show()

Circle is fixed at North Pole, it won't move

swatchai · Accepted Answer · 2019-07-12T17:16:56.447

1

The coordinates of the circle's center must be the projection coordinates. So, a coordinate transformation is required.

Here is the relevant code:

# Note: lat = 72, lon = 100
# proj = ccrs.NorthPolarStereo(central_longitude=lon)

projx1, projy1 = proj.transform_point(lon, lat, ccrs.Geodetic()) #get proj coord of (lon,lat)
ax.add_patch(mpatches.Circle(xy=[projx1, projy1], radius=r_ortho, color='red', \
                            alpha=0.3, transform=proj, zorder=30))

The output plot will be:

Alternate Solution

Since the projection in use is conformal, a Tissot Indicatrix plot on it will always be a circle. So that, the circle you need can be represented with an Indicatrix. Here is the code you can try:

ax.tissot(rad_km=r_ortho/1000, lons=lon, lats=lat, n_samples=48, color='green', \
      alpha=0.3, zorder=31)

Edit 1

Code to replace the errorneous function:

import pyproj

def compute_radius(val_degree):
    """
    Compute surface distance in meters for a given angular value in degrees
    """
    geod84 = pyproj.Geod(ellps='WGS84')
    lat0, lon0 = 0, 90
    _, _, dist_m = geod84.inv(lon0, lat0,  lon0+val_degree, lat0)
    return dist_m

compute_radius(1)  # get: 111319.49079327357 (meters)

edited Jul 12 '19 at 17:16

answered Jul 12 '19 at 08:08

swatchai

17,400
3
39
58

thanks for the solution. However, I would expect the radius `r` to be expressed in degrees. If I decrease `r=1` for example, I still get a huge circle. Am I missing something here? – ouranos Jul 12 '19 at 10:40
Indeed, the transformation of the radius into projection coordinates doesn't work as I would expect. Could you please shed some light here? – ouranos Jul 12 '19 at 10:54
@ouranos The extents of the plot relate to the value of `pad_radius` which depends on `r`. When you change the value of `r` from 20 to 1 the new plot will covers smaller earth extent (but equal plot size). While the circle size remains the same on the map, but relatively, it becomes smaller. – swatchai Jul 12 '19 at 11:17
that is what I expected to see. But it doesn't happen so. Setting `r=1` I would expect a radius of 1 degree length, but the circle still cover almost the entire Siberia. The difference in radius length from 20 to 1 doesn't make sense to me yet. – ouranos Jul 12 '19 at 14:32
If it works normally for you, could you please paste the entire code in your answer? I may be missing something to be able to reproduce it. – ouranos Jul 12 '19 at 14:33
1

Function `compute_radius()` does not give correct values as expected. – swatchai Jul 12 '19 at 16:37
1

@ouranos I added the code for `compute_radius()` that should work for you. – swatchai Jul 12 '19 at 17:22

Drawing circles with Cartopy in NorthPolarStereo projection

1 Answers1

Linked