Python matplotlib/cartopy data on curvilinear rotated pole grid#

Software requirements:

  • Python 3

  • matplotlib

  • cartopy

  • numpy

  • xarray

Example script#

import xarray as xr
import numpy as np
import os, sys
import matplotlib.pyplot as plt
import cartopy
import as ccrs

def read_data(file_name):
    Read netcdf file and return arrays: lat, lon, var, pole_longitude, pole_latitude.
    ds = xr.open_dataset(file_name)
    var  = ds.tas[0,:,:]
    rlat = ds.rlat[:]
    rlon = ds.rlon[:]
    pole = ds.rotated_pole

        if hasattr(pole,'grid_north_pole_longitude'):
            px = pole.attrs['grid_north_pole_longitude']
        if hasattr(pole,'grid_north_pole_latitude'):
            py = pole.attrs['grid_north_pole_latitude']
        print('Unexpected error:', sys.exc_info()[0])

    return rlon, rlat, var, px, py

def main():
    Main function.
    colormap = "RdYlBu_r"

    fname = ''

    # read file content and return relevant variables
    rlon, rlat, var, pole_lon, pole_lat = read_data(fname)

    ax = plt.axes(projection=ccrs.PlateCarree())

    ax.set_extent([-46, 70, 20, 75], crs=ccrs.PlateCarree())

    ax.add_feature(cartopy.feature.OCEAN, color='white', zorder=0)
    ax.add_feature(cartopy.feature.LAND, color='lightgray',zorder=0,
                   linewidth=0.5, edgecolor='black')

    ax.gridlines(draw_labels=True, linewidth=0.5, color='gray',
                 xlocs=range(-180,180,15), ylocs=range(-90,90,15))

    ax.coastlines(resolution='50m', linewidth=0.3, color='black')

    ax.set_title('Py: de-rotated grid (cartopy)', fontsize=10, fontweight='bold')

    crs = ccrs.RotatedPole(pole_longitude=pole_lon, pole_latitude=pole_lat)
    ax.contourf(rlon, rlat, var, levels=15, cmap=colormap, transform=crs)

    plt.savefig('plot_Python_curvilinear_grid_derotated_1.png', bbox_inches='tight', dpi=200)

if __name__ == '__main__':