CODEWITHSUNDEEP
projectioncartopy
ocean_phd
ocean_phd

Reputation: 23

Cartopy figure for high latitude with edges parallel to latitude and longitude, e.g., not rectangular

I'm trying to create a Cartopy map for the sub-polar region around Iceland. What I would like is a non-rectangular figure where the edges are parallel to the lines of longitude and latitude, like this figure created using PyGMT:

Map between 45 and 70 degrees north and 45 and 5 degrees west with edges of map parallel to longitude and latitude

I've tried various Cartopy projections, but all result in a rectangular figure, e.g.,

import matplotlib.pyplot as plt
import cartopy
import cartopy.crs as ccrs

fig = plt.figure()
proj = ccrs.LambertConformal(central_longitude=-25, central_latitude=58.0)
ax = plt.axes(projection = proj)

ax.set_extent((-45, -5, 45, 70))
ax.gridlines()
ax.add_feature(cartopy.feature.LAND, zorder=1, edgecolor='black')

There are reasons for not using PyGMT (I want to plot surface velocities using quiver, plus the extensive learning curve), so I'd like to know if it's possible to achieve the same result in cartopy.

Thanks

Upvotes: 1

Views: 1443

Answers (2)

swatchai
swatchai

Reputation: 18812

The answer by @Rutger_Kassies is great. However there is an alternative that the readers should consider if he/she wants to try a different approach.

import cartopy.crs as ccrs
import cartopy
import matplotlib.pyplot as plt
import matplotlib.path as mpath

# The lat-long proj
noProj = ccrs.PlateCarree()
# The projection of the map
myProj = ccrs.LambertConformal(central_longitude=-25, central_latitude=58.0)
myProj._threshold = myProj._threshold/20.  # Set for higher precision of the projection

ax = plt.axes(projection=myProj)

# This plots parallel and meridian arcs around a target area that will be used ...
#  as the map boundary
[ax_hdl] = ax.plot([-45, -5, -5, -45, -45], [45, 45, 70, 70, 45],
         color='black', linewidth=0.5, marker='none',
         transform=noProj)
# Get the `Path` of the plot
tx_path = ax_hdl._get_transformed_path()
path_in_data_coords, _ = tx_path.get_transformed_path_and_affine()

# Use the path's vertices to create a polygon
polygon = mpath.Path( path_in_data_coords.vertices )
ax.set_boundary(polygon) #This masks-out unwanted part of the plot

ax.gridlines(draw_labels=True, x_inline=False, y_inline=False)
ax.add_feature(cartopy.feature.OCEAN, linewidth=.3, color='lightblue')
ax.add_feature(cartopy.feature.LAND, zorder=1, edgecolor='black')
ax.title.set_text("Meridians and Parallels as Boundary")
plt.show()

merpar

You can change some parameters in the code, for example, the type of arcs that are used as the map's boundary.

The second plot is obtained by changing these parts of code:

1. `transform=noProj`  to  
   `transform=ccrs.Geodetic()`
2. `ax.title.set_text("Meridians and Parallels as Boundary")` to
   `ax.title.set_text("Great-circle Arcs as Boundary")`

mer_gcc

I believe that when the top edge of the map reaches high latitude, using parallel of latitude as the boundary is not optimum. Straight line may be better, but in some situation great circle arc should be considered.

Upvotes: 1

Rutger Kassies
Rutger Kassies

Reputation: 64463

You can use the set_boundary method of an axes for this. When specifying it as lon/lat, for a different projection, you should sample some points accross the boundary to approximate the true curvature of the projection (compared to lon/lat). The example below takes 20 points on each edge.

Note that the shape of this boundary can be anything you want, it doesn't have to match the projection or lon/lat lines etc.

import matplotlib.pyplot as plt
import matplotlib.path as mpath
import cartopy
import cartopy.crs as ccrs
import numpy as np

proj = ccrs.LambertConformal(central_longitude=-25, central_latitude=58.0)

fig, axs = plt.subplots(
    1,2, figsize=(8, 3), facecolor="w", 
    subplot_kw=dict(projection=proj),
)

n = 20
aoi = mpath.Path(
    list(zip(np.linspace(-45,-5, n), np.full(n, 70))) + \
    list(zip(np.full(n, -5), np.linspace(70, 45, n))) + \
    list(zip(np.linspace(-5, -45, n), np.full(n, 45))) + \
    list(zip(np.full(n, -45), np.linspace(45, 70, n)))
)

axs[1].set_boundary(aoi, transform=ccrs.PlateCarree())
    
for ax in axs:
    ax.set_extent((-45, -5, 45, 70))
    ax.add_feature(cartopy.feature.LAND, zorder=1, edgecolor='k')
    gl = ax.gridlines(
        draw_labels=True, rotate_labels=False,
        x_inline=False, y_inline=False,
    )

enter image description here

Upvotes: 1

Related Questions