Lesson 28: CartoPy Basics

This lesson is modified from Introduction to Cartopy by Project Pythia, and Maps with CartoPy by Research Computing in Earth Science

---

Overview

CartoPy is a Python library that provides cartographic tools for creating maps and visualizing geospatial data with Matplotlib. By then end of this lesson you will be able to

1. Explain core Cartopy concepts such as map projections and GeoAxes

2. . Creating regional maps with raster data (where each cell contains a single value) or vector data (represting geographic features as points, lines, and polygons)

3. Create a map figure with multiple subplots

---

Installation and import

We import the main libraries of Cartopy: Coordinate Reference Systems (crs) and Feature.

• CRS defines how coordinates are projected onto a flat surface.

• Features in Cartopy represent elements such as coastlines, borders, and rivers that can be added to maps for better visualization and context.

#pip install cartopy

from cartopy import crs as ccrs
from cartopy import feature as cfeature
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import xarray as xr

#  Suppress warnings issued by Cartopy when downloading data files
import warnings
warnings.filterwarnings('ignore')

---

1. Basic concepts

Python integratez with libraries and tools for data analysis, visualization, and machine learning. This allows you to create end-to-end workflows. An important step in the workflow is map plotting.

1.1 Extend Matplotlib’s axes into georeferenced GeoAxes

In previous tutorials, we learned that a figure in Matplotlib consists of a Figure object and a list of one or more Axes objects (subplots). By importing cartopy.crs, we gain access to Cartopy’s Coordinate Reference System, which offers various geographical projections. Specifying one of these projections for an Axes object transforms it into a GeoAxes object, effectively georeferencing the subplot.

1.2 Create a map with a specified projection

Here, we will generate a GeoAxes object with the PlateCarree projection. The PlateCarree projection is a global lat-lon map projection where each point is evenly spaced in terms of degrees. The term “Plate Carré” is French for “flat square.”

# Create a figure with a specified size
fig = plt.figure(figsize=(11, 8.5))

# Add a subplot with PlateCarree projection
ax = plt.subplot(1, 1, 1, projection=ccrs.PlateCarree(central_longitude=-75))
ax.set_title("A Geo-referenced subplot, Plate Carree projection")

Text(0.5, 1.0, 'A Geo-referenced subplot, Plate Carree projection')

Although the figure seems empty, it has been georeferenced using a map projection using Cartopy’s crs (coordinate reference system) class. We can now add in cartographic features to our subplot in the form of shapefiles. One such cartographic feature is coastlines, which can be added to our subplot using the callable GeoAxes method simply called coastlines.

# Method 1 for pre-defined features of Natural Earht datasets
ax.coastlines(resolution='110m', color='black')
fig

Info

To get the figure to display again with the features that we've added since the original display, just type the name of the Figure object in its own cell.

1.3 Add cartographic features to the map

Cartopy provides other cartographic features via its features class, which was imported at the beginning of this page, under the name cfeature. These cartographic features are represented as data stored in shapefiles. When these cartographic features are called for the first time in a script or notebook, the shapefiles are downloaded from Natural Earth Data. A list of the different Natural Earth shapefiles can be accessed from CartoPy documentation including a list of CartoPy pre-defined features (e.g. .BORDERS and .STATES).

To incorporate these features into our subplot, we can utilize the add_feature method, which permits the specification of attributes (e.g., line style) through arguments (e.g. linestyle=dashed), similar to Matplotlib’s plot method.

Let us include borders and U.S. state lines in our subplot.

# Adding borders and states to the plot
ax.add_feature(cfeature.BORDERS, linewidth=0.5, edgecolor='black')
ax.add_feature(cfeature.STATES, linewidth=0.3, edgecolor='brown')

fig

In the above example, we used predefined features like cfeature.BORDERS, but we can use any Natural Earth dataset by creating a CartoPy instance with cfeature.NaturalEarthFeature() as explained further below.

1.4 Add a background image

CartoPy has only one pre-created background image that you can use with .stock_img method. You can add other background images such as etopo toppgraphy or NASA Blue marble using background_img() method. In addition to background images, you can use ax.imshow() in CartoPy to display any image on your map as detailed below.

Let us add the pre-created background image stock_img:

#Add the only pre-created background image in CartoPy
ax.stock_img()

fig

1.5 Cartopy’s map projections

The crs class of CartoPy offers a variety of map projections for visualizing geospatial data such as

1. Plate Carree that represents the globe as a rectangular grid with equal spacing in all directions, and thus exaggerates the spatial extent of regions closer to the poles.

2. Mercator, which is a cylindrical map projection that preserves straight lines, making it useful for navigation purposes.

3. Robinson that balances distortions in size and shape, suitable for world maps.

4. Orthographic, a perspective projection where the globe is projected onto a tangent plane, creating a realistic view of the Earth.

5. Mollweide, an equal-area projection that distorts shape and direction while preserving area, making it ideal for thematic maps and often used with global satellite mosaics

In the Plate Carree (Equirectangular) projection, the meridians (lines of longitude) and parallels (lines of latitude) are straight and intersect at right angles similar to Mecator. This characteristic contributes to the simplicity of the Plate Carree projection, as the grid lines form a regular rectangular grid on the map. However. the distortion increases as you move towards the poles, making the higher latitudes appear larger than they actually are.

projections = [ccrs.PlateCarree(),
               ccrs.Mercator(),
               ccrs.Robinson(),
               ccrs.Orthographic(),
               ccrs.Mollweide(),
              ]

for proj in projections:
    plt.figure()
    ax = plt.axes(projection=proj)
    ax.stock_img()
    ax.coastlines()
    ax.add_feature(cfeature.BORDERS, linewidth=0.5, edgecolor='blue');
    ax.set_title(f'{type(proj)}');

2. Create regional maps

For this example, let us plot a map for the Gulf of Mexico. We will :

1. set extend to the Gulf of Mexico lat and lon coordinates

2. add features such as coastlines, boarder, states, land, rivers, ocean, and bathymetry

3. plot features such as markers for major cities, a track line, and study area boundaries

4. plot 2d raster data such as sea surface height (zos)

5. Add and georeference an image to our map

Let us do this part by part.

2.1 Cartopy’s set_extent method

Using Plate Carree projection, we will use Cartopy’s set_extent method to restrict the map coverage to a Florida view.

Here we set the domain, which defines the geographical region to be plotted that is used in the set_extent call. Since these coordinates are expressed in degrees, they correspond to a PlateCarree projection.

#Select plot extend: Lon_West, Lon_East, Lat_South, Lat_North in Degrees
extent = [-75, -100, 18, 31] 
cLon = (extent[0] + extent[1]) / 2
cLat = (extent[2] + extent[3]) / 2

#Select projection
projPC = ccrs.PlateCarree()

Source:https://www.satsig.net/lat_long.htm

Let us do this regional plot including plotting the latitude and longitude lines, and adding a title.

# Create figure and axes
fig = plt.figure(figsize=(11, 8.5))
ax = plt.subplot(1, 1, 1, projection=projPC)

# Setting the extent of the plot
ax.set_extent(extent, crs=projPC)

# Add plot title and gridlines
ax.set_title('Plate Carree')
ax.gridlines(draw_labels=True, linewidth=1, color='gray', alpha=0.5, linestyle='--')

<cartopy.mpl.gridliner.Gridliner at 0x2420f344290>

Info

The call to the subplot method use PateCarree projection, the calls to set_extent use PlateCarree. This ensures that the values we passed into set_extent will be transformed from degrees into the values appropriate for the projection we use for the map. Here it makes no difference since the map projection is also PateCarree.

2.2 Adding features from Natural Earth Datasets

In CartoPy, you can add features from Natural Earth datasets to a map in various ways.

# Method 1 for Cartopy pre-defined features of Natural Earht datasets
ax.add_feature(cfeature.COASTLINE.with_scale('10m'),color='black')`,

# Method 2 beyond Cartopy pre-defined features of Natural Earht datasets
ax.add_feature(cfeature.NaturalEarthFeature(
    category='physical', 
    name='coastline', 
    scale='10m', 
    facecolor='none', 
    edgecolor='black'))

The first method is suitable for CartoPy’s pre-defined features including BORDERS, COASTLINE, LAKES, LAND, OCEAN, RIVERS, and STATES, while the second method allows you to incorporate additional shapefiles from Natural Earth beyond the pre-defined features. When using pre-defined features, you can specify resolutions of ‘10m’, ‘50m’, or ‘110m. For custom Natural Earth features, ensure to verify the available resolutions for the feature you intend to add by referring to the Natural Earth features website.

2.2.1 Add some predefined features

For illustration, let us use the first method to add some pre-defined features such as boarders, land, rivers, and lakes. Some cartographical features are predefined as constants in the cartopy.feature package. The resolution of these features depends on the amount of geographical area in your map, specified by set_extent.

# Set background color to water color
ax.set_facecolor(cfeature.COLORS['water'])

# Add Land with zorder=0 meaning that feature order will be backward
ax.add_feature(cfeature.LAND.with_scale('10m'), zorder=20) 

# Add Coastlines, and Lakes above Land features
ax.add_feature(cfeature.LAKES.with_scale('10m'), color='aqua', zorder=21) 
ax.add_feature(cfeature.COASTLINE.with_scale('10m'), color='gray', zorder=21) 

# Add Boarders feature on the top
ax.add_feature(cfeature.BORDERS, 
               edgecolor='red',facecolor='none', linewidth=2, zorder=22) 

fig

Note:

For high-resolution Natural Earth shapefiles such as this, while we could add Cartopy's OCEAN feature, it currently takes much longer to render on the plot. Instead, we take the strategy of first setting the facecolor of the entire subplot to match that of water bodies in Cartopy. When we then layer on the LAND feature, pixels that are not part of the LAND shapefile remain in the background as water facecolor, which is the same color as the OCEAN.

By setting different zorder values for different map features, you control their drawing order on the plot canvas. Features with higher zorder values are drawn on top of those with lower values. The above code adds land features and sets the zorder to 20, meaning that these land features will be drawn below features with higher zorder values such as lakes and coastlines with zorder of 21 to ensure that lakes and coastlines are drawn above land features. The code adds borders to the plot with zorder of 22, placing the borders on top of all other features. We will later assign zorder values of 10 to 19 to ocean features to ensure that they drawn below land featurs.

2.2.2 Add shapefiles from Natural Earth

Instead of using pre-defined features, here we will use shapefiles from Natural Earth. This requires much more code than previous examples. We first need to create new objects associated shapefiles using the NaturalEarthFeature method, which is part of the Cartopy feature class. Then we can use add_feature to add the new objects to our new map.

The general syntax is

ax.add_feature(cfeature.NaturalEarthFeature(
    category='physical', 
    name='coastline', 
    scale='10m', 
    facecolor='none', 
    edgecolor='black'))

The categorycan be:

• ‘physical’ with name such as ‘coastline’, ‘land’, ‘rivers_lake_centerlines’, etc.

• ‘cultural’ with name such as ‘state’, ‘urban polygons’, admin_1_states_provinces, etc.

and the scale can be 10, 50, or 110.

To find the feature category, name and scale:

• select the feature that you want from Natural Earth Data

• and obtain the category, name and scale from the feature download link.

For example, for rivers in North America, the download link is:

http//www.naturalearthdata.com/download/10m/physical/ne_10m_rivers_north_america.zip

Accordingly, we can set category='physical', name='rivers_north_america', and scale='10m', and so on.

Let us add rivers and supplementary North America rivers. In addition, let us add few bathymetry lines using a loop.

# Add rivers
ax.add_feature(cfeature.NaturalEarthFeature(
    'physical',
    'rivers_lake_centerlines',
    scale='10m',
    edgecolor='blue',
    facecolor='none')
               ,zorder=21)

# Add additional for North America
ax.add_feature(cfeature.NaturalEarthFeature(
    'physical',
    'rivers_north_america',
    scale='10m',
    edgecolor='blue',
    facecolor='none',
    alpha=0.5),
               zorder=21)

# Plot 6 bathymetry lines using a loop
# Names of bathymetry line 
bathymetry_names = ['bathymetry_L_0', 'bathymetry_K_200', 
                    'bathymetry_J_1000', 'bathymetry_I_2000', 
                    'bathymetry_H_3000', 'bathymetry_G_4000'] 

# Very light blue to light blue generated by ChatGPT 3.5 Turbo
colors = [(0.9, 0.95, 1), (0.7, 0.85, 1), 
          (0.5, 0.75, 1), (0.3, 0.65, 1), 
          (0.1, 0.55, 1), (0, 0.45, 1), ]  

# Loop to plot each bathymetry line
for bathymetry_name, color in zip(bathymetry_names, colors):
    ax.add_feature(
        cfeature.NaturalEarthFeature(
            category='physical', 
            name=bathymetry_name, 
            scale='10m', 
            edgecolor='none', 
            facecolor=color), 
        zorder = 11)

fig

Warning

Be patient; when plotting a small geographical area, the high-resolution "10m" shapefiles are used by default. As a result, these plots take longer to create, especially if the shapefiles are not yet downloaded from Natural Earth. Similar issues can occur whenever a `GeoAxes` object is transformed from one coordinate system to another.

2.3 plot user-defined features

Let us as study area boundaries, a track line, and markers for major cities.

# Plot site area on GeoAxes given lon and lat coordinates
ax.plot([-85.0, -85.0], [26.5, 28.0], color='red', linewidth=1, zorder = 15)  
ax.plot([-85.0, -82.85], [28.0, 28.0], color='red', linewidth=1, zorder = 15)
ax.plot([-85.0, -81.9], [26.5, 26.5], color='red', linewidth=1, zorder = 15)


# Plot Track 026 line on GeoAxex and annotate with text using axes cooridnates
# Deinfe start and end points of the line
points = np.array([[-85.7940, 28.4487], [-83.6190, 23.8411]]) 

# Plot the line
ax.plot(points[:,0],points[:,1], color='black', linewidth=2, zorder = 15)

# Annotate the line
ax.text(0.6, 0.6,  # Coordinates are in axes coordinate rather than data coordinates  
        'Track 026',  # Specify text content
        transform=ax.transAxes,  # Set coordinate system for coordinates, which is the coordinate system of the axes itself 
        color='black',  # Set color of the text
        fontsize=10,  # Set font size of the text
        ha='center',  # Set horizontal alignment of the text
        va='center',  # Set vertical alignment of the text
        rotation=300, # Rotate text
        zorder = 15,  # Plot order of the feature
       )

# Add names and locations of major cities
# Define markers for each location as a dictionary (city coordinates generated by ChatGPT 3.5 Turbo)
locations_coordinates = {
    "Tampa": [27.9506, -82.4572],
    "Sarasota": [27.3364, -82.5307],
    "Fort Myers": [26.6406, -81.8723],
    "Naples": [26.142, -81.7948]}

# Use a loop to add markers and text labels slightly away from each other for each location 
for location, coordinates in locations_coordinates.items():
    #Plot marker
    ax.plot(coordinates[1],  # lon coordinate of locatin in Degrees
            coordinates[0],  # lat coordinate of locatin in Degrees
            transform=ccrs.PlateCarree(),  # Set coordinate system for coordinates, which is Plate Carree
            marker='o',      # Marker style for the plot
            color='black',     # Color of the markers
            markersize=3,    # Size of the markers   
            zorder = 22,  # Plot order of the feature
           )
    #Annote each city
    ax.text(coordinates[1] + 0.1,   # lon coordinate of locatin in Degrees (slightly adjusted away from marker)
            coordinates[0] + 0.1,   # lat coordinate of locatin in Degrees (slightly adjusted away from marker)
            location,               # Text to display (city name)
            transform=ccrs.PlateCarree(),   # Specify coordinate system for coordinates, which is Plate Carree
            color='black',           # Color of the text
            fontsize=8,              # Font size of the text
            zorder = 22,  # Plot order of the feature
           )

fig

Note:

The transform=ax.transAxes parameter is used to specify the coordinate system of the provided coordinates. In this case, ax.transAxes refers to the coordinate system of the axes itself. When you use ax.transAxes, the coordinates (0,0) correspond to the lower left corner of the axes and (1,1) correspond to the upper right corner of the axes, regardless of the data limits or scales within the plot. By specifying transform=ax.transAxes, you are indicating that the coordinates (0.6, 0.6) are relative to the size of the axes, ensuring that the text 'Track 026' is positioned at 60% of the width and 60% of the height of the axes, regardless of the data being plotted.

2.4 Plot 2d raster data

Plotting data on a Cartesian grid is like using the PlateCarree projection, where meridians and parallels are straight lines with constant spacing. We can plot our raster data as a contour map. Similar to what we did before, we use the transform keyword argument in the contour-plotting method to specify the current projection of our raster data, which will be transformed into the projection of our subplot method. In this particular case, we do not need to use use the transform keyword since the project of our raster data is Plate Carree, which is the same as our subplot projection.

We can read raster data from a NetCDF file using xr.open_dataset(), or generate our data. Let us generate some data that we can add to our figure.

# Create coordinates for the four labels
dates = pd.date_range('2024-01-03', periods=2)  # Generating a time range for the data
depths = [10, 20, 30]  # List of depth values
latitudes = np.linspace(24.5,30.5,80)  # latitude values (y-values)
longitudes =np.linspace(-87,-77,100)  # longitude values (x-values)

# Create a 4D Numpy Array representing values  for a 2x3x4x5 grid
size=(len(dates), len(depths), len(latitudes), len(longitudes))
zos = np.round(np.random.uniform(-1, 1, size= size), 2) #sea surface height [m]
tos = np.round(np.random.uniform(26, 30, size=size), 2) #sea surface temperature [C]
sos = np.round(np.random.uniform(33, 36, size=size), 2) #sea surface salinty (PSU) 
                 
# Create a DataArray with zos values
ds = xr.Dataset({'zos': (('date', 'depth', 'lat', 'lon'), zos),
                'tos': (('date', 'depth', 'lat', 'lon'), tos),
                'sos': (('date', 'depth', 'lat', 'lon'), sos),},
    coords={'date': dates, 'depth': depths, 'lat': latitudes, 'lon': longitudes})

# Add metadata to the DataArray
ds.zos.attrs['Description'] = "Sea surface height above the geoid (m)"
ds.tos.attrs['Description'] = "Sea surface temperature (°C)"
ds.sos.attrs['Description'] = "Sea surface salinty (PSU)" 

# Display ds
ds

<xarray.Dataset> Size: 1MB
Dimensions:  (date: 2, depth: 3, lat: 80, lon: 100)
Coordinates:
  * date     (date) datetime64[ns] 16B 2024-01-03 2024-01-04
  * depth    (depth) int64 24B 10 20 30
  * lat      (lat) float64 640B 24.5 24.58 24.65 24.73 ... 30.35 30.42 30.5
  * lon      (lon) float64 800B -87.0 -86.9 -86.8 -86.7 ... -77.2 -77.1 -77.0
Data variables:
    zos      (date, depth, lat, lon) float64 384kB -0.29 0.98 ... -0.67 -1.0
    tos      (date, depth, lat, lon) float64 384kB 28.37 27.69 ... 28.62 27.14
    sos      (date, depth, lat, lon) float64 384kB 35.79 35.92 ... 34.88 35.74

xarray.Dataset

• Dimensions:
  • date: 2
  • depth: 3
  • lat: 80
  • lon: 100
• Coordinates: (4)
  • date
    (date)
    datetime64[ns]
    2024-01-03 2024-01-04

    array(['2024-01-03T00:00:00.000000000', '2024-01-04T00:00:00.000000000'],
          dtype='datetime64[ns]')

  • depth
    (depth)
    int64
    10 20 30

    array([10, 20, 30])

  • lat
    (lat)
    float64
    24.5 24.58 24.65 ... 30.42 30.5

    array([24.5     , 24.575949, 24.651899, 24.727848, 24.803797, 24.879747,
           24.955696, 25.031646, 25.107595, 25.183544, 25.259494, 25.335443,
           25.411392, 25.487342, 25.563291, 25.639241, 25.71519 , 25.791139,
           25.867089, 25.943038, 26.018987, 26.094937, 26.170886, 26.246835,
           26.322785, 26.398734, 26.474684, 26.550633, 26.626582, 26.702532,
           26.778481, 26.85443 , 26.93038 , 27.006329, 27.082278, 27.158228,
           27.234177, 27.310127, 27.386076, 27.462025, 27.537975, 27.613924,
           27.689873, 27.765823, 27.841772, 27.917722, 27.993671, 28.06962 ,
           28.14557 , 28.221519, 28.297468, 28.373418, 28.449367, 28.525316,
           28.601266, 28.677215, 28.753165, 28.829114, 28.905063, 28.981013,
           29.056962, 29.132911, 29.208861, 29.28481 , 29.360759, 29.436709,
           29.512658, 29.588608, 29.664557, 29.740506, 29.816456, 29.892405,
           29.968354, 30.044304, 30.120253, 30.196203, 30.272152, 30.348101,
           30.424051, 30.5     ])

  • lon
    (lon)
    float64
    -87.0 -86.9 -86.8 ... -77.1 -77.0

    array([-87.      , -86.89899 , -86.79798 , -86.69697 , -86.59596 , -86.494949,
           -86.393939, -86.292929, -86.191919, -86.090909, -85.989899, -85.888889,
           -85.787879, -85.686869, -85.585859, -85.484848, -85.383838, -85.282828,
           -85.181818, -85.080808, -84.979798, -84.878788, -84.777778, -84.676768,
           -84.575758, -84.474747, -84.373737, -84.272727, -84.171717, -84.070707,
           -83.969697, -83.868687, -83.767677, -83.666667, -83.565657, -83.464646,
           -83.363636, -83.262626, -83.161616, -83.060606, -82.959596, -82.858586,
           -82.757576, -82.656566, -82.555556, -82.454545, -82.353535, -82.252525,
           -82.151515, -82.050505, -81.949495, -81.848485, -81.747475, -81.646465,
           -81.545455, -81.444444, -81.343434, -81.242424, -81.141414, -81.040404,
           -80.939394, -80.838384, -80.737374, -80.636364, -80.535354, -80.434343,
           -80.333333, -80.232323, -80.131313, -80.030303, -79.929293, -79.828283,
           -79.727273, -79.626263, -79.525253, -79.424242, -79.323232, -79.222222,
           -79.121212, -79.020202, -78.919192, -78.818182, -78.717172, -78.616162,
           -78.515152, -78.414141, -78.313131, -78.212121, -78.111111, -78.010101,
           -77.909091, -77.808081, -77.707071, -77.606061, -77.505051, -77.40404 ,
           -77.30303 , -77.20202 , -77.10101 , -77.      ])

• Data variables: (3)
  • zos
    (date, depth, lat, lon)
    float64
    -0.29 0.98 0.66 ... 0.34 -0.67 -1.0

    Description :

    Sea surface height above the geoid (m)

    array([[[[-0.29,  0.98,  0.66, ...,  0.53,  0.47, -0.12],
             [-0.16,  0.32,  0.2 , ..., -0.74,  0.44,  0.03],
             [-0.63,  0.57,  0.58, ...,  0.75,  0.56, -0.35],
             ...,
             [-0.48, -0.27, -0.24, ...,  0.96, -0.44,  0.71],
             [ 0.88, -0.75, -0.26, ..., -0.52, -0.26, -0.66],
             [ 0.87,  0.86, -0.25, ...,  1.  , -0.11, -0.54]],
    
            [[ 0.26, -0.66,  0.98, ..., -0.23,  0.67,  0.8 ],
             [-0.07,  0.06, -0.39, ..., -0.91,  0.81, -0.8 ],
             [-0.4 , -0.37, -0.59, ...,  0.15,  0.78, -0.56],
             ...,
             [ 0.43,  0.75, -0.92, ...,  0.74, -0.62, -0.69],
             [ 0.88, -0.39,  0.24, ..., -0.96,  0.39, -0.03],
             [-0.75,  0.95,  0.51, ...,  0.23,  0.77, -0.84]],
    
            [[ 0.07,  0.65, -0.83, ..., -0.19,  0.2 , -0.99],
             [ 0.05, -0.18, -0.14, ...,  0.14,  0.85, -0.99],
             [-0.54,  0.94,  0.2 , ..., -0.16, -0.46, -0.93],
             ...,
    ...
             ...,
             [-0.47,  0.89, -0.37, ...,  0.75, -0.57,  0.28],
             [-0.97, -0.24, -0.88, ...,  0.82,  0.61, -0.96],
             [ 0.33, -0.63, -0.14, ...,  0.34,  0.41,  0.52]],
    
            [[-0.12,  0.03, -0.43, ...,  0.56,  0.21, -0.19],
             [-0.79,  0.55, -0.15, ..., -0.82, -0.73, -0.37],
             [-0.28,  0.89,  0.74, ..., -0.41,  0.46, -0.78],
             ...,
             [ 0.52, -0.15,  0.08, ...,  0.54, -0.97, -0.88],
             [-0.89,  0.78, -0.52, ..., -0.77, -0.56,  0.49],
             [ 0.66,  0.6 ,  0.9 , ..., -0.51,  0.25, -0.17]],
    
            [[-0.05, -0.01,  0.43, ...,  0.81,  0.45,  0.35],
             [-0.38,  0.42,  0.46, ...,  0.41, -0.22, -0.31],
             [-0.42,  0.6 , -0.8 , ..., -0.74,  0.91,  0.28],
             ...,
             [-0.47,  0.64,  0.96, ..., -0.06, -0.15,  0.78],
             [ 0.42, -0.07,  0.92, ..., -0.25,  0.8 , -0.15],
             [ 0.65, -0.08, -0.03, ...,  0.34, -0.67, -1.  ]]]])

  • tos
    (date, depth, lat, lon)
    float64
    28.37 27.69 27.81 ... 28.62 27.14

    Description :

    Sea surface temperature (°C)

    array([[[[28.37, 27.69, 27.81, ..., 27.15, 27.44, 28.18],
             [27.39, 27.01, 26.67, ..., 29.75, 29.42, 26.08],
             [26.55, 29.89, 29.39, ..., 28.13, 27.41, 29.13],
             ...,
             [27.75, 27.23, 29.08, ..., 28.6 , 28.79, 28.44],
             [27.34, 28.63, 29.66, ..., 28.28, 29.78, 26.89],
             [28.71, 27.19, 28.53, ..., 27.01, 29.99, 27.56]],
    
            [[29.75, 29.87, 26.78, ..., 26.28, 29.81, 26.95],
             [26.61, 28.52, 26.57, ..., 29.1 , 28.64, 29.3 ],
             [26.5 , 27.22, 27.65, ..., 29.86, 29.81, 27.73],
             ...,
             [26.17, 29.02, 26.14, ..., 29.94, 26.24, 29.31],
             [28.88, 27.23, 26.25, ..., 29.89, 27.08, 28.36],
             [29.45, 28.7 , 26.24, ..., 29.33, 26.02, 27.75]],
    
            [[27.9 , 27.3 , 29.14, ..., 29.53, 26.48, 29.59],
             [27.96, 27.94, 29.52, ..., 27.43, 26.22, 27.89],
             [26.92, 29.83, 26.06, ..., 28.86, 28.57, 27.48],
             ...,
    ...
             ...,
             [27.6 , 28.92, 26.01, ..., 27.21, 29.29, 26.22],
             [26.26, 26.7 , 26.8 , ..., 26.16, 27.55, 27.35],
             [29.78, 26.79, 27.96, ..., 29.15, 27.62, 27.82]],
    
            [[28.33, 29.47, 27.68, ..., 27.59, 27.38, 28.61],
             [26.25, 27.7 , 29.58, ..., 27.35, 29.36, 28.86],
             [26.49, 28.11, 29.67, ..., 27.81, 28.98, 26.33],
             ...,
             [29.86, 27.49, 28.99, ..., 27.69, 28.57, 29.43],
             [26.05, 29.66, 26.73, ..., 29.09, 26.06, 26.69],
             [29.88, 26.33, 27.8 , ..., 29.95, 28.79, 27.9 ]],
    
            [[29.12, 26.13, 27.53, ..., 27.56, 29.02, 26.68],
             [26.72, 28.4 , 29.85, ..., 26.91, 29.66, 29.32],
             [26.01, 28.16, 27.83, ..., 29.88, 28.18, 28.47],
             ...,
             [29.34, 26.19, 29.37, ..., 27.54, 26.61, 27.93],
             [29.95, 28.03, 26.91, ..., 28.14, 26.97, 27.89],
             [27.94, 27.4 , 27.58, ..., 28.19, 28.62, 27.14]]]])

  • sos
    (date, depth, lat, lon)
    float64
    35.79 35.92 33.62 ... 34.88 35.74

    Description :

    Sea surface salinty (PSU)

    array([[[[35.79, 35.92, 33.62, ..., 33.84, 33.82, 35.51],
             [34.5 , 33.76, 33.98, ..., 35.6 , 34.87, 34.38],
             [33.04, 33.7 , 35.9 , ..., 35.26, 35.29, 34.51],
             ...,
             [35.84, 34.81, 35.56, ..., 34.89, 35.59, 35.43],
             [33.42, 33.05, 35.36, ..., 33.99, 34.93, 35.05],
             [33.52, 33.06, 35.86, ..., 34.74, 33.94, 33.56]],
    
            [[35.1 , 34.09, 35.65, ..., 34.96, 33.93, 34.47],
             [33.47, 34.48, 34.17, ..., 35.52, 34.76, 33.08],
             [35.04, 35.59, 34.13, ..., 35.75, 35.07, 34.49],
             ...,
             [34.39, 33.92, 33.82, ..., 35.9 , 35.12, 34.01],
             [35.34, 33.38, 33.14, ..., 35.67, 34.27, 34.17],
             [33.89, 33.71, 35.19, ..., 33.27, 34.67, 34.51]],
    
            [[35.  , 35.21, 34.12, ..., 35.29, 34.85, 33.69],
             [35.37, 33.08, 34.41, ..., 35.6 , 34.07, 33.15],
             [35.3 , 35.05, 33.13, ..., 34.24, 33.51, 35.93],
             ...,
    ...
             ...,
             [34.33, 33.47, 33.12, ..., 35.77, 33.97, 33.28],
             [34.4 , 35.17, 34.54, ..., 35.32, 35.6 , 34.36],
             [34.31, 33.4 , 34.28, ..., 33.92, 35.81, 33.96]],
    
            [[34.28, 35.97, 33.08, ..., 33.82, 33.02, 33.55],
             [35.49, 33.58, 35.06, ..., 33.91, 35.27, 35.78],
             [33.56, 33.22, 34.35, ..., 35.08, 35.76, 34.61],
             ...,
             [34.1 , 35.87, 35.22, ..., 34.77, 33.38, 35.41],
             [33.88, 35.37, 33.68, ..., 34.56, 33.59, 35.4 ],
             [33.62, 34.43, 34.63, ..., 34.33, 33.28, 34.44]],
    
            [[34.92, 35.77, 35.87, ..., 33.14, 34.94, 35.85],
             [34.56, 34.  , 35.91, ..., 34.28, 34.62, 34.14],
             [33.11, 34.73, 34.41, ..., 34.24, 34.  , 34.  ],
             ...,
             [35.96, 33.37, 33.09, ..., 33.83, 35.26, 34.6 ],
             [33.72, 34.43, 33.86, ..., 33.83, 35.1 , 33.69],
             [33.4 , 35.37, 34.48, ..., 34.68, 34.88, 35.74]]]])

• Indexes: (4)
  • date
    PandasIndex

    PandasIndex(DatetimeIndex(['2024-01-03', '2024-01-04'], dtype='datetime64[ns]', name='date', freq='D'))

  • depth
    PandasIndex

    PandasIndex(Index([10, 20, 30], dtype='int64', name='depth'))

  • lat
    PandasIndex

    PandasIndex(Index([              24.5, 24.575949367088608, 24.651898734177216,
           24.727848101265824, 24.803797468354432, 24.879746835443036,
           24.955696202531644, 25.031645569620252,  25.10759493670886,
            25.18354430379747, 25.259493670886076, 25.335443037974684,
           25.411392405063292,   25.4873417721519, 25.563291139240505,
           25.639240506329113,  25.71518987341772,  25.79113924050633,
           25.867088607594937, 25.943037974683545, 26.018987341772153,
            26.09493670886076,  26.17088607594937, 26.246835443037973,
            26.32278481012658,  26.39873417721519, 26.474683544303797,
           26.550632911392405, 26.626582278481013,  26.70253164556962,
            26.77848101265823, 26.854430379746837, 26.930379746835442,
            27.00632911392405, 27.082278481012658, 27.158227848101266,
           27.234177215189874, 27.310126582278482,  27.38607594936709,
           27.462025316455698, 27.537974683544306,  27.61392405063291,
           27.689873417721518, 27.765822784810126, 27.841772151898734,
           27.917721518987342,  27.99367088607595, 28.069620253164558,
           28.145569620253163, 28.221518987341774,  28.29746835443038,
           28.373417721518987, 28.449367088607595, 28.525316455696203,
            28.60126582278481,  28.67721518987342, 28.753164556962027,
            28.82911392405063, 28.905063291139243, 28.981012658227847,
           29.056962025316455, 29.132911392405063,  29.20886075949367,
            29.28481012658228, 29.360759493670887, 29.436708860759495,
             29.5126582278481,  29.58860759493671, 29.664556962025316,
           29.740506329113924,  29.81645569620253,  29.89240506329114,
           29.968354430379748, 30.044303797468356, 30.120253164556964,
           30.196202531645568,  30.27215189873418, 30.348101265822784,
           30.424050632911392,               30.5],
          dtype='float64', name='lat'))

  • lon
    PandasIndex

    PandasIndex(Index([             -87.0,  -86.8989898989899, -86.79797979797979,
            -86.6969696969697,  -86.5959595959596,  -86.4949494949495,
           -86.39393939393939, -86.29292929292929,  -86.1919191919192,
            -86.0909090909091, -85.98989898989899, -85.88888888888889,
           -85.78787878787878, -85.68686868686869, -85.58585858585859,
           -85.48484848484848, -85.38383838383838, -85.28282828282828,
           -85.18181818181819, -85.08080808080808, -84.97979797979798,
           -84.87878787878788, -84.77777777777777, -84.67676767676768,
           -84.57575757575758, -84.47474747474747, -84.37373737373737,
           -84.27272727272727, -84.17171717171718, -84.07070707070707,
           -83.96969696969697, -83.86868686868686, -83.76767676767676,
           -83.66666666666667, -83.56565656565657, -83.46464646464646,
           -83.36363636363636, -83.26262626262626, -83.16161616161617,
           -83.06060606060606, -82.95959595959596, -82.85858585858585,
           -82.75757575757575, -82.65656565656566, -82.55555555555556,
           -82.45454545454545, -82.35353535353535, -82.25252525252526,
           -82.15151515151516, -82.05050505050505, -81.94949494949495,
           -81.84848484848484, -81.74747474747474, -81.64646464646465,
           -81.54545454545455, -81.44444444444444, -81.34343434343434,
           -81.24242424242425, -81.14141414141415, -81.04040404040404,
           -80.93939393939394, -80.83838383838383, -80.73737373737374,
           -80.63636363636364, -80.53535353535354, -80.43434343434343,
           -80.33333333333333, -80.23232323232324, -80.13131313131314,
           -80.03030303030303, -79.92929292929293, -79.82828282828282,
           -79.72727272727273, -79.62626262626263, -79.52525252525253,
           -79.42424242424242, -79.32323232323232, -79.22222222222223,
           -79.12121212121212, -79.02020202020202, -78.91919191919192,
           -78.81818181818181, -78.71717171717172, -78.61616161616162,
           -78.51515151515152, -78.41414141414141, -78.31313131313131,
           -78.21212121212122, -78.11111111111111, -78.01010101010101,
            -77.9090909090909,  -77.8080808080808, -77.70707070707071,
           -77.60606060606061,  -77.5050505050505,  -77.4040404040404,
            -77.3030303030303, -77.20202020202021,  -77.1010101010101,
                        -77.0],
          dtype='float64', name='lon'))

• Attributes: (0)

Let us add a slice of this zos data to our map, add a colorbar for zos, and add a label to the color bar.

# Select the data for the specific time and depth
data_slice = ds.zos.isel(date=0, depth=0)  

dataplot = ax.contourf(data_slice.lon,   # X-coordinates for the contour plot
                       data_slice.lat,   # Y-coordinates for the contour plot
                       data_slice,       # Data values to plot
                       cmap='coolwarm',  # Colormap for the plot
                       zorder=12,         # Drawing order for the plot above bathymetry and below land and study area
                       transform=ccrs.PlateCarree()  #Not needed as both our data and figure are Plate Carree
                      )

# Add a colorbar
colorbar = plt.colorbar(dataplot,  # Create a colorbar for dataplot data
                        orientation='horizontal',  # Set the orientation of the colorbar to horizontal
                        pad=0.06,  # Set the padding between the colorbar and the plot
                        shrink=0.8)  # Adjust the size of the colorbar to 80% of the default size

# Add a label to the colorbar with a specified font size
colorbar.set_label(ds.zos.attrs['Description'],  fontsize=12)


fig

<Figure size 640x480 with 0 Axes>

2.5 Adding and georeferencing an image to a map

We can add a satellite image on the map if we know the coordinates of the image edges.

# Read image 
img = plt.imread('Data/Track026.png')
img_extent = (-90, -81, 23, 31) #lon E to W and Lat S to N

# Add image to map 
ax.imshow(img,              # Image to plot
          origin='upper',   # Origin is upper left corner for png image
          extent=img_extent, 
          transform=ccrs.PlateCarree(), # Specify coordinate system for image coordinates 
          zorder=12,        # Drawing order for image is above raster data but below land
         )

fig

---------------------------------------------------------------------------
FileNotFoundError                         Traceback (most recent call last)
Cell In[15], line 2
      1 # Read image 
----> 2 img = plt.imread('Data/Track026.png')
      3 img_extent = (-90, -81, 23, 31) #lon E to W and Lat S to N
      5 # Add image to map 

File ~\AppData\Local\miniconda3\Lib\site-packages\matplotlib\pyplot.py:2607, in imread(fname, format)
   2603 @_copy_docstring_and_deprecators(matplotlib.image.imread)
   2604 def imread(
   2605         fname: str | pathlib.Path | BinaryIO, format: str | None = None
   2606 ) -> np.ndarray:
-> 2607     return matplotlib.image.imread(fname, format)

File ~\AppData\Local\miniconda3\Lib\site-packages\matplotlib\image.py:1512, in imread(fname, format)
   1505 if isinstance(fname, str) and len(parse.urlparse(fname).scheme) > 1:
   1506     # Pillow doesn't handle URLs directly.
   1507     raise ValueError(
   1508         "Please open the URL for reading and pass the "
   1509         "result to Pillow, e.g. with "
   1510         "``np.array(PIL.Image.open(urllib.request.urlopen(url)))``."
   1511         )
-> 1512 with img_open(fname) as image:
   1513     return (_pil_png_to_float_array(image)
   1514             if isinstance(image, PIL.PngImagePlugin.PngImageFile) else
   1515             pil_to_array(image))

File ~\AppData\Local\miniconda3\Lib\site-packages\PIL\ImageFile.py:132, in ImageFile.__init__(self, fp, filename)
    128 self.decodermaxblock = MAXBLOCK
    130 if is_path(fp):
    131     # filename
--> 132     self.fp = open(fp, "rb")
    133     self.filename = os.fspath(fp)
    134     self._exclusive_fp = True

FileNotFoundError: [Errno 2] No such file or directory: 'Data/Track026.png'

3. A figure with multiple subplots

Let put everything together and create a figure with three subplots for the three variables and one with an image, using a loop.

# 1) Select plot extend: Lon_West, Lon_East, Lat_South, Lat_North in Degrees
extent = [-75, -100, 18, 31] 
cLon = (extent[0] + extent[1]) / 2
cLat = (extent[2] + extent[3]) / 2

# 2) Select projection
projPC = ccrs.PlateCarree()

# 3) Create figure
fig = plt.figure(figsize=(11, 8.5))
plt.subplots_adjust(hspace=0.1, wspace=0.2)  # Adjust the height and width spaces between subplots

# 4) Loop to plot variables 
variables = ['zos','tos','sos', 'image']
titles = ['(a) ', '(b)', '(c)', '(d)']
colormaps = ['coolwarm', 'jet', 'viridis', '']

for index, (variable, title, colormap) in enumerate(zip(variables, titles, colormaps)):

    # 5) Create Axes     
    ax = plt.subplot(2, 2, index+1, projection=projPC)
    
    # 6) Setting the extent of the plot
    ax.set_extent(extent, crs=projPC)
    
    # 7) Add plot title 
    ax.set_title(title, fontsize=10, loc='left')
    
    # 8) Add gridlines (we need to cutomize grid lines to improve visualization)                              
    gl = ax.gridlines(draw_labels=True, xlocs=np.arange(extent[1],extent[0], 5), ylocs= np.arange(extent[2], extent[3], 5),
                      linewidth=1, color='gray', alpha=0.5, linestyle='--',  zorder=40)
    gl.right_labels = False  # 
    gl.top_labels = False    # Set the font size for x-axis labels
    gl.xlabel_style = {'size': 10}  # Set the font size for x-axis labels
    gl.ylabel_style = {'size': 10}  # Set the font size for y-axis labels
    
    # 9) Adding Natural Earth features
    ## Pre-defined features
    ax.set_facecolor(cfeature.COLORS['water'])  #Blackgroundwater color for the ocean
    ax.add_feature(cfeature.LAND.with_scale('10m'), zorder=20) 
    ax.add_feature(cfeature.LAKES.with_scale('10m'), color='aqua', zorder=21) 
    ax.add_feature(cfeature.COASTLINE.with_scale('10m'), color='gray', zorder=21) 
    ax.add_feature(cfeature.BORDERS, edgecolor='red',facecolor='none', linewidth=2, zorder=22) 

    ## General features (not pre-defined)
    ax.add_feature(cfeature.NaturalEarthFeature(
        'physical', 'rivers_lake_centerlines', scale='10m', edgecolor='blue', facecolor='none') ,zorder=21)
    ax.add_feature(cfeature.NaturalEarthFeature(
        'physical','rivers_north_america', scale='10m', edgecolor='blue', facecolor='none', alpha=0.5),zorder=21)

    ## Multiple features with a loop
    bathymetry_names = ['bathymetry_L_0', 'bathymetry_K_200', 'bathymetry_J_1000', 
                        'bathymetry_I_2000',  'bathymetry_H_3000', 'bathymetry_G_4000'] 
    colors = [(0.9, 0.95, 1), (0.7, 0.85, 1),  (0.5, 0.75, 1), (0.3, 0.65, 1),  (0.1, 0.55, 1), (0, 0.45, 1), ]  
    for bathymetry_name, color in zip(bathymetry_names, colors):
        ax.add_feature(cfeature.NaturalEarthFeature(
            category='physical', name=bathymetry_name, scale='10m', edgecolor='none', facecolor=color), zorder = 11)

    # 10) Adding user-defined features 
    ## Plot site area on GeoAxes given lon and lat coordinates in first subplot 
    if index == 0:
        ax.plot([-85.0, -85.0], [26.5, 28.0], color='black', linewidth=1, zorder = 15)  
        ax.plot([-85.0, -82.85], [28.0, 28.0], color='black', linewidth=1, zorder = 15)
        ax.plot([-85.0, -81.9], [26.5, 26.5], color='black', linewidth=1, zorder = 15)

    ## Plot Track 026 line on GeoAxex and annotate with text using axes cooridnates in second subplot
    if index == 1: 
        points = np.array([[-85.7940, 28.4487], [-83.6190, 23.8411]]) 
        ax.plot(points[:,0],points[:,1], color='black', linewidth=2, zorder = 15)
        ax.text(0.58, 0.6, 'Track 026', transform=ax.transAxes, 
                color='black', fontsize=10,  ha='center',  va='center',  rotation=300, zorder = 15)

    ## Adding names and locations of major cities
    locations_coordinates = {"Tampa": [27.9506, -82.4572], 
                             #"Sarasota": [27.3364, -82.5307], 
                             #"Fort Myers": [26.6406, -81.8723], 
                             "Naples": [26.142, -81.7948]}
    for location, coordinates in locations_coordinates.items():
        ax.plot(coordinates[1], coordinates[0],  
                transform=ccrs.PlateCarree(),  marker='o', color='black', markersize=3, zorder = 22)
        ax.text(coordinates[1] + 0.1, coordinates[0] + 0.1, 
                location, transform=ccrs.PlateCarree(), color='black', fontsize=8, zorder = 22)

    # 11) Adding raster data and colorbar for the first three subplots
    if index < 3: 
        data_slice = ds[variable].isel(date=0, depth=0)  
        dataplot = ax.contourf(data_slice.lon, data_slice.lat, data_slice, 
                               cmap=colormap, zorder=12, transform=ccrs.PlateCarree())
        colorbar = plt.colorbar(dataplot, orientation='horizontal', pad=0.1, shrink=1)  
        colorbar.set_label(ds[variable].attrs['Description'],  fontsize=10)

    # 12) Adding an image for the fourth subplot 
    else:
        img = plt.imread('Data/Track026.png')
        img_extent = (-90, -81, 23, 31) #lon E to W and Lat S to N
        imgplot = ax.imshow(img, origin='upper', extent=img_extent, transform=ccrs.PlateCarree(), zorder=12)
        
        # adding invisible colorbar for the layout of the subplots remains aligned visually
        # you can also do this by adjusting subplot size or using GridSpec as suggested by ChatGPT 3.5 Turbo
        cbar = plt.colorbar(imgplot, orientation='horizontal', pad=0.1, shrink=1)  
        cbar.ax.set_visible(False)

Creating publication-quality figures with subplots is time-consuming, as it involves balancing automation with attention to detail. Factors contributing to this include precise layout adjustments, customization of plot elements, ensuring consistency, and handling complex visualizations. This iterative process of refinement is important for achieving a polished final figure. Libraries like Seaborn or Plotly can streamline complex visualization creation with less manual effort.

4. Saving a figure

CartoPy is built on top of Matplotlib. When saving a figure in Matplotlib, you can specify the format using the savefig() function. For example, to save a figure as a PNG file, you can use:

plt.savefig('figure.png', format='png')

Commonly used image formats in Matplotlib include:

• PNG for sharp edges and solid colors with lossless compression

• TIFF for high-quality images with lossless compression

• JPEG for smooth gradients or photographs with small file sizes

• SVG for scalable plots without quality loss as a vector format

• PDF for embedding in documents

You can choose the image format based on your specific requirements such as quality, scalability, and intended use of the plot.

Let’s save our last image with subplots.

#Save the figure as an image file
plt.savefig('output_image.png', dpi=300, bbox_inches='tight')

<Figure size 640x480 with 0 Axes>

---

Summary

• Cartopy allows for the georeferencing of Matplotlib Axes objects.

• Cartopy’s crs class supports a variety of map projections.

• Cartopy’s feature class allows for a variety of cartographic features to be overlaid on a georeferenced plot or subplot.

---

Resources and References

1. Cartopy Documentation

2. Full list of projections in Cartopy

3. Maps with Cartopy (Ryan Abernathey)

4. Map Projections (GeoCAT)

5. NCAR xdev Cartopy Tutorial Video

6. Python Tutorial Seminar Series - Cartopy