Python matplotlib example: Overlay regional on global grid

Software requirements:

  • Python 3

  • Numpy

  • matplotlib

  • cartopy

  • xarray

Projection: - Orthographic

Grids:

  1. The first grid is a coarse global grid (CMIP5)

    • grid type: lonlat

    • points=18432 (192x96)

    • lon : 0 to 358.125 by 1.875 degrees_east  circular

    • lat : -88.57217 to 88.57217 degrees_north

  2. The second grid is a high resolution regional grid (CORDEX EUR-11)

    • grid type: curvilinear

    • points=174688 (424x412)

    • lon : -44.59386 to 64.96438 degrees_east

    • lat : 21.98783 to 72.585 degrees_north

Script:

#!/usr/bin/env python
# coding: utf-8
#
# DKRZ matplotlib example: Draw two grids with different resolutions
#
# Used packages:
#     - numpy
#     - xarray
#     - matplotlib
#     - cartopy
#
# Projection:
#     - Orthographic
#
# Grids:
#
# 1. The first grid is a coarse global grid (CMIP5)
#
#     grid type: lonlat
#     points=18432 (192x96)
#     lon : 0 to 358.125 by 1.875 degrees_east  circular
#     lat : -88.57217 to 88.57217 degrees_north
#
# 2. The second grid is a high resolution regional grid (CORDEX EUR-11)
#
#     grid type: curvilinear
#     points=174688 (424x412)
#     lon : -44.59386 to 64.96438 degrees_east
#     lat : 21.98783 to 72.585 degrees_north
#

import numpy as np
import xarray as xr

import cartopy.crs as ccrs
import cartopy.feature as cfeature
import cartopy.util as cutil

import matplotlib.pyplot as plt
import matplotlib.patches as mpatches

def main():

    #-- Read global data:

    fname1 = "orog_fx_MPI-ESM-LR_rcp26_r0i0p0.nc"    #-- orography
    mname1 = "sftlf_fx_MPI-ESM-LR_rcp26_r0i0p0.nc"   #-- land area fraction

    #--  open file and read variables
    f1    =  xr.open_dataset(fname1)
    var1  =  f1.variables["orog"][:,:]

    mask1 =  xr.open_dataset(mname1)
    lsm1  =  mask1.variables["sftlf"][:,:]
    lsm1  =  np.where(lsm1 > 0.5, lsm1, 0)

    land_only1 = np.where(lsm1 > 0.5, var1, -101)    #-- min contour value - 1

    lat   =  f1.variables["lat"][:]
    lon   =  f1.variables["lon"][:]
    dlat  =  lat[1]-lat[0]
    dlon  =  lon[1]-lon[0]

    #-- Add cyclic points:

    cyclic_data, cyclic_lon = cutil.add_cyclic_point(land_only1, coord=lon)

    #-- Read regional data:

    fname2 = "orog_EUR-11_MPI-M-MPI-ESM-LR_historical_r0i0p0_CLMcom-CCLM4-8-17_v1_fx.nc"
    mname2 = "sftlf_EUR-11_MPI-M-MPI-ESM-LR_historical_r0i0p0_CLMcom-CCLM4-8-17_v1_fx.nc"

    f2    = xr.open_dataset(fname2)
    var2  = f2.variables["orog"][:,:]

    mask2 = xr.open_dataset(mname2)
    lsm2  = mask2.variables["sftlf"][:,:]
    lsm2  = np.where(lsm2 > 0.5, lsm2, 0)

    land_only2 = np.where(lsm2 > 0.5, var2, -101)    #-- min contour value - 1

    lat2d = f2.variables["lat"][:,:]
    lon2d = f2.variables["lon"][:,:]

    nlat  = len(lat2d[:,0])
    nlon  = len(lon2d[0,:])

    #-- Define edges of regional data:

    lon_val_lower = lon2d[0,:]
    lon_val_right = lon2d[:,nlon-1]
    lon_val_left  = lon2d[:,0]
    lon_val_upper = lon2d[nlat-1,:]

    lat_val_lower = lat2d[0,:]
    lat_val_right = lat2d[:,nlon-1]
    lat_val_left  = lat2d[:,0]
    lat_val_upper = lat2d[nlat-1,:]

    #-- Generate the data for the edges of the regional grid

    line_lons = np.append([lon_val_upper], [lon_val_right[::-1]])
    line_lons = np.append([line_lons], [lon_val_lower])
    line_lons = np.append([line_lons], [lon_val_left])
    line_lats = np.append([lat_val_upper], [lat_val_right[::-1]])
    line_lats = np.append([line_lats], [lat_val_lower])
    line_lats = np.append([line_lats], [lat_val_left])

    polyline = np.column_stack([line_lons, line_lats])

    #-- Create the color mesh plot:

    projection=ccrs.Orthographic(central_latitude=50.0, central_longitude=10.0)

    fig, ax = plt.subplots(subplot_kw=dict(projection=projection), figsize=(10,9))

    ax.set_global()

    #-- add coastal outlines
    ax.add_feature(cfeature.COASTLINE.with_scale('50m'), linewidth=0.5)

    gl = ax.gridlines(draw_labels=True, linewidth=0.5, color='k', zorder=3)
    gl.xlabel_style = {'size':10}
    gl.ylabel_style = {'size':10}
    gl.top_labels   = True
    gl.right_labels = True

    #-- plot the title string
    plt.title('Overlay: regional on global grid (orography)')

    #-- define color map
    cmap = 'terrain'

    #-- create the color mesh plot
    edgecolor = [0.0, 0.0, 0.0, 0.2] #'k'
    linewidth = 0.005

    cnf1  = ax.pcolormesh(cyclic_lon, lat, cyclic_data,
                          cmap=cmap,
                          vmin=-800,
                          vmax=3000,
                          edgecolor=edgecolor,
                          linewidth=linewidth,
                          transform=ccrs.PlateCarree())

    cnf2 = ax.pcolormesh(lon2d, lat2d, land_only2,
                         cmap=cmap,
                         vmin=-800,
                         vmax=3000,
                         transform=ccrs.PlateCarree())

    #-- add a polyline around the regional grid
    lw, ec, fc = 1, 'k', 'y'     #-- linewidth, edgecolor, facecolor

    ax.add_patch(mpatches.Polygon(polyline,
                                  closed=False,
                                  fill=False,
                                  linewidth=lw,
                                  edgecolor=ec,
                                  facecolor=fc,
                                  transform=ccrs.Geodetic()))

    #-- add a color bar
    cbar_ax = fig.add_axes([0.94, 0.25, 0.015, 0.6], autoscalex_on=True)  #-- x,y,w,h
    cbar    = fig.colorbar(cnf1, cax=cbar_ax, orientation='vertical')
    plt.setp(cbar.ax.get_xticklabels()[::2], visible=False)
    cbar.set_label('[m]')

    #-- save the plot in PNG format
    plt.savefig('grid_resolutions_overlay.png', bbox_inches='tight', dpi=100)

if __name__ == '__main__':
    main()

Result:

image0