Search code examples
pythonmatplotlibcartopy

Cartopy fails with small regional plots in polar regions

I want to make simple static maps for use in journal papers. I work in the Arctic and I make a lot of map images that show equipment layout, vessel tracks, and source and receiver locations. I can't do it in CARTOPY. For example, 1 deg of latitude (74 to 75N) and 2.5 degrees of longitude (-92.5 to -90.0) at a mid-latitude of say 74.5N. You can't get the coastline to work properly. The map is often empty, but it should show a portion of the coastline of Devon Island, NU. If I make the plot a bigger region (something like 30deg by 30deg), it works, but you will see that the coordinates displayed in the graph window don't line up properly. The X, Y values match the graph axes, but the lat, lon values are in parentheses are shifted. In the worst case, the lat, lons come out as 0.00n deg or even nearly half a world away.

I've tried multiple ways of invoking the extent. Different projections. Nothing seems to work.

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
from cartopy.feature import NaturalEarthFeature
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
import numpy as np
import matplotlib.ticker as mticker


# the limits of the map
# extent = (-100., -50.0, 60.0, 80.0)  # try this, you'll get a shifted plot on northern Alaska
extent = (-92.5, -90.0, 74.0, 75.0)   # try this, you'll get a blank plot. axes shifted badly.

# set the projection type
c_lon, c_lat = (extent[0] + extent[1])/2., (extent[2] + extent[3])/2.
proj = ccrs.PlateCarree(central_longitude=c_lon)
# proj = ccrs.Mercator(c_lon)    # I've tried these as well
# proj = ccrs.Orthographic(c_lon, c_lat)

gax = plt.axes(projection=proj)
gax.set_extent(extent)
gax.set_ylim((extent[2], extent[3]))
gax.set_xlim((extent[0], extent[1]))

# now add the coastline. This only works for big maps. Small regions fail.
coastline_10m = NaturalEarthFeature(category='physical', name='coastline', \
                    facecolor='none', scale='10m')
gax.add_feature(coastline_10m, edgecolor='gray')

# draw a grid with labelled lat and lon. Suppress ticklabels on the top and right.
gl = gax.gridlines(crs=proj, draw_labels=True) # only works with PlateCarree()
gl.xlabels_top = None
gl.ylabels_right = False

# now we put labels on the X and Y axes. You have to move these around manually.
gax.text(-0.2, 0.55, 'Latitude [Deg]', va='bottom', ha='center',
        rotation='vertical', rotation_mode='anchor',
        transform=gax.transAxes)
gax.text(0.5, -0.12, 'Longitude [Deg]', va='bottom', ha='center',
        rotation='horizontal', rotation_mode='anchor',
        transform=gax.transAxes)

# approximately correct for the aspect ratio
plt.gca().set_aspect(1.0/(np.cos(np.pi*(extent[2] + extent[3])/(2.*180.))))
plt.show()

macOS Catalina, Anaconda python3.8.3, IPython 7.19.0, cartopy 0.17 (the highest version supported. Anaconda says 0.18, but it installs 0.17).

Solution

  • Some obvious errors:

    1. set_extent() needs option crs=ccrs.PlateCarree()
    2. set_xlim() and set_ylim need data (map projection) coordinates.

    And set_xlim() and set_ylim change what you have done with set_extent() in the previous line. They must be used correctly. Most of your cases, they should not be used.