Python use AlbersEqualArea projection#
In this example the use of the Albers Equal Area map projection is demonstrated. The data used here is the CMIP6 MPI-ESM1-2-LR ssp585 sea ice concentration.
Software requirements:
Python 3
matplotlib
cartopy
python-cdo
Example script#
siconc_AlbersEqualArea_projection.py
#!/usr/bin/env python
# coding: utf-8
'''
DKRZ example
Polar plot of sea ice concentration data using Albers Equal Area projection and
a satellite image as background
-------------------------------------------------------------------------------
Description
In this example we demonstrate how to draw the CMIP6 sea ice concentration data
on top of the satellite image of the Earth.
The map projection is Albers Equal Area, see
https://scitools.org.uk/cartopy/docs/latest/reference/projections.html#albersequalarea
-------------------------------------------------------------------------------
Data
CMIP6 MPI-ESM1-2-LR SSP5-8.5 sea ice concentration. <br>
Further info URL : <br>
https://furtherinfo.es-doc.org/CMIP6.MPI-M.MPI-ESM1-2-LR.ssp585.none.r10i1p1f1
-------------------------------------------------------------------------------
Satellite image of the Earth
The image of the Earth - The Blue Marble - with a high resolution of 0.1 degrees
(3600x1800 pixels) for August 2004 can be downloaded from the NASA Earth
Observatory download page <br>
https://neo.gsfc.nasa.gov/view.php?datasetId=BlueMarbleNG&date=2004-08-01
To display the satellite image as map background the `CARTOPY_USER_BACKGROUNDS`
environment variable which points to the **images.json** file has to be set.
The satellite image has to be stored in the same directory as **images.json**
that will be read by Cartopy's `cartopy.mpl.geoaxes.GeoAxes.background_image()`
method.
-------------------------------------------------------------------------------
2024 copyright DKRZ licensed under CC BY-NC-SA 4.0
(https://creativecommons.org/licenses/by-nc-sa/4.0/deed.en)
-------------------------------------------------------------------------------
'''
import os
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
import matplotlib.path as mpath
import matplotlib.colors as mcolors
import matplotlib.ticker as mticker
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
from cdo import Cdo
cdo = Cdo()
def main():
#------------------------------------------------------------------------------
# Function: transform_to_proj()
#
# Transform function to prevent the line wrapping of data on curvilinear grids
# when `contourf` or `contour` plot style is used.
#------------------------------------------------------------------------------
def transform_to_proj(longitude, latitude,
projection=ccrs.NorthPolarStereo(),
transform=ccrs.PlateCarree()):
lat = latitude.values.ravel()
lon = longitude.values.ravel()
X, Y = [], []
for itr in range(len(lat)):
x, y = projection.transform_point(lon[itr], lat[itr], transform)
X.append(x)
Y.append(y)
shapes = latitude.shape
X = np.reshape(X, (shapes[0], shapes[1]))
Y = np.reshape(Y, (shapes[0], shapes[1]))
return X, Y
#------------------------------------------------------------------------------
# Path to background images and `images.json` file
#
# Tell Cartopy where to find the `images.json` file that keeps the map background
# image details.
#
# See https://scitools.org.uk/cartopy/docs/latest/reference/generated/cartopy.mpl.geoaxes.GeoAxes.html#cartopy.mpl.geoaxes.GeoAxes.background_img
#------------------------------------------------------------------------------
os.environ['CARTOPY_USER_BACKGROUNDS'] = os.environ['HOME']+'/Python/Backgrounds'
#------------------------------------------------------------------------------
# Sea ice concentration data
#
# Read the sea ice concentration data but use only the data on the northern hemisphere.
# This makes the computations faster. To extract the data and compute the seasonal
# means we use the Python bindings of CDO, python-cdo.
#
# Compute the seasonal means but for DJF left out the first and the last one,
# because of incompleteness. Set seasmean values equal zero to missing value.
#------------------------------------------------------------------------------
input_file = os.environ['HOME']+'/data/CMIP6/LR/SeaIce/Runs/ssp585/LR_SImon_EM_SICONC.nc'
input_cmd = f'-seasmean ' + \
f'-seltimestep,3/-2 ' + \
f'-sellonlatbox,0.,360.,30.,90. {input_file}'
ds_season = cdo.setctomiss('0.', input=input_cmd, returnXDataset=True)
variable = ds_season.siconc
scen_name = 'SSP5-8.5'
season_dict = {'mam': 0, 'jja': 1, 'son': 2, 'djf': 3}
# Define the order of the seasons in a dictionary and extract the season data.
season_dict = dict(mam=0, jja=1, son=2, djf=3)
var_mam = variable.isel(time=range(season_dict['mam'], ds_season.time.size, 4))
var_jja = variable.isel(time=range(season_dict['jja'], ds_season.time.size, 4))
var_son = variable.isel(time=range(season_dict['son'], ds_season.time.size, 4))
var_djf = variable.isel(time=range(season_dict['djf'], ds_season.time.size, 4))
season_data = [var_mam, var_jja, var_son, var_djf]
# Choose the time step and extract the variable data we want to show. Set values
# equal zero to NaN.
use_season = 'MAM'
use_timestep = 0
var = season_data[use_timestep].isel(time=use_timestep)
year = var.time.dt.year.data
#------------------------------------------------------------------------------
# Plotting
#
# Choose the plot style whether to use `contour`, `contourf`, or `pcolormesh`.
# Default plotstyle is `pcolormesh`.
#------------------------------------------------------------------------------
allowed_plotstyles = ['pcolormesh', 'contour', 'contourf']
use_plotstyle = 'contourf'
# Set the data coordinate reference system and the map projection data. Only
# the northern part for latitude 70° to 88° and longitude -90° to 90° should be
# displayed.
data_crs = ccrs.PlateCarree()
projection = ccrs.AlbersEqualArea(central_longitude=0., central_latitude=79.,
standard_parallels=(70, 88))
# To prevent the line wrapping for data on a curvilinear grid using contour or
# contourf plot style, we have to compute the data transformation to the
# NorthPolarStereo projection. This takes some time.
if use_plotstyle in ('contour', 'contourf'):
X, Y = transform_to_proj(ds_season.longitude, ds_season.latitude,
projection=projection,
transform=data_crs)
# Colormap and normalization for values 0 to 100% in 5% steps.
bounds = np.arange(0, 105, 5)
norm = mcolors.BoundaryNorm(boundaries=bounds, ncolors=256)
# Create the Albers Equal Area plot with the satellite image in the background.
plt.switch_backend('agg')
fig, ax = plt.subplots(figsize=(12,9), subplot_kw=dict(projection=projection))
lon_min, lon_max = -90., 90.
lat_min, lat_max = 60., 90.
ax.set_extent([lon_min, lon_max, lat_min, lat_max], crs=data_crs)
ax.background_img(name='BlueMarbleBright', resolution='high')
ax.add_feature(cfeature.BORDERS.with_scale('50m'), lw=1., ls='--', ec='k', zorder=3)
ax.add_feature(cfeature.COASTLINE.with_scale('50m'), lw=0.5, ec='k', zorder=4)
#-- check plot style
if not use_plotstyle in allowed_plotstyles:
print(f'Warning: use_plotstyle={use_plotstyle} is not allowed. Set use_plotstyle to pcolormesh')
use_plotstyle = 'pcolormesh'
#-- set common plot settings
plot_kw = dict(cmap='Blues_r', norm=norm)
if use_plotstyle == 'pcolormesh':
plot = ax.pcolormesh(var.longitude, var.latitude, var.data, transform=data_crs, **plot_kw)
else: #-- Note: in this case do not use transform keyword!
plot = getattr(ax, use_plotstyle)(X, Y, var.data, levels=20, **plot_kw)
#-- add a colorbar
cbar = plt.colorbar(plot, shrink=0.6, pad=0.04, orientation='horizontal',
drawedges=True, aspect=50, anchor=(0.5, 1.), ticks=bounds)
cbar.set_ticklabels(bounds)
cbar.set_label('%', labelpad=0)
#-- cutting out the semi-circular area
vertices = [(longitude, lat_min) for longitude in np.arange(lon_min , lon_max, 1)] + \
[(longitude, lat_max) for longitude in np.arange(lon_max, lon_min-1, -1)]
boundary = mpath.Path(vertices)
ax.set_boundary(boundary, transform=data_crs)
#-- add grid lines and axes tick labels
gl = ax.gridlines(draw_labels=True, color='k', ls='--', lw=0.2)
#-- x-tick labels
gl.xlocator = mticker.FixedLocator(np.arange(-90, 90, 10))
gl.xformatter = LONGITUDE_FORMATTER
#-- y-tick labels
gl.xlabel_style = dict(weight='normal', color='black', rotation=0,)
gl.ylocator = mticker.FixedLocator(np.arange(50, 90, 10))
gl.yformatter = LATITUDE_FORMATTER
gl.ylabel_style = dict(weight='normal', color='black', rotation=0,)
#-- add title and other annotations
dy = 0.03
tx = fig.text(0.15, 0.82+dy, 'CMIP6', fontsize=16, weight='normal')
tx = fig.text(0.51, 0.82+dy, 'Sea Ice Concentration', fontsize=16, ha='center', weight='normal')
tx = fig.text(0.51, 0.78+dy, scen_name+'/ MPI-ESM1-2-LR', fontsize=10, weight='normal', ha='center')
tx = fig.text(0.82, 0.83, f'{year}', fontsize=16, weight='normal', ha='center')
tx = fig.text(0.82, 0.81, f'{use_season}', fontsize=11, weight='normal', ha='center')
tx = fig.text(0.75, 0.24, '© 2024 CC BY-NC-SA 4.0 4.0: DKRZ', fontsize=7)
#-- save the figure to PNG
fig.savefig('plot_siconc_'+scen_name+'_AlbersEqualArea_proj.png',
bbox_inches='tight', facecolor='white', dpi=100)
if __name__ == '__main__':
main()
Plot result#
