NCL_lcnative_1.py

This script illustrates the following concepts:

• Drawing contours over a map using a native lat,lon grid

• Drawing filled contours over a Lambert Conformal map

• Zooming in on a particular area on a Lambert Conformal map

• Subsetting a color map

• Using all three Cartopy Lambert projections to find the best fit for visualization

• Implementing best practices for choosing contour color scheme. The original NCL version uses a rainbow color map which can create problems for individuals impacted by color blindness and is not black and white printer friendly. To overcome this, we used the “Blues_r” color scheme from Matplotlib. More information on best practices can be found here.

See following URLs to see the reproduced NCL plot & script:

• Original NCL script: https://www.ncl.ucar.edu/Applications/Scripts/lcnative_1.ncl

• Original NCL plot: https://www.ncl.ucar.edu/Applications/Images/lcnative_1_lg.png

Import packages:

import numpy as np
import xarray as xr
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker

import geocat.datafiles as gdf

Read in data:

# Open a netCDF data file using xarray default engine and load the data into xarrays
ds = xr.open_dataset(gdf.get("netcdf_files/pre.8912.mon.nc"),
                     decode_times=False)

# Extract a slice of the data
t = ds.pre.isel(time=0)

Plot:

plt.figure(figsize=(14, 14))


def Plot(row, col, pos, proj, title):
    '''
    Args:

        row (:class: 'int'):
            number of rows necessary for subplotting of visualizations
        col (:class: 'int'):
            number of columns necessary for subplotting
        pos (:class: 'int'):
            position of visualization in m x n subplot
        proj (:class: 'cartopy.crs'):
            which projection to visualize
        title (:class: 'str'):
            center title of respective visualization
    '''
    # Generate axes using Cartopy and draw coastlines
    ax = plt.subplot(row, col, pos, projection=proj)
    ax.set_extent((28, 57, 20, 47), crs=ccrs.PlateCarree())
    ax.coastlines(linewidth=0.5)

    # Use ax "gridlines" function to draw lat/lon markers on projections
    gl = ax.gridlines(draw_labels=True,
                      dms=False,
                      x_inline=False,
                      y_inline=False)
    gl.top_labels = True
    gl.right_labels = True
    gl.xlines = False
    gl.ylines = False
    gl.xlocator = mticker.FixedLocator([30, 35, 40, 45, 50, 55])
    gl.ylocator = mticker.FixedLocator([20, 25, 30, 35, 40, 45])
    gl.xlabel_style = {'rotation': 0}
    gl.ylabel_style = {'rotation': 0}
    '''
    When using certain types of projections in Cartopy, you may find that there
    is not a direct 1-to-1 projection similarity. When looking at the three Lambert
    projections offered, you will notice the closest match to the NCL projection
    is actually the Lambert Cylindrical projection. This is due to NCL having certain
    "smoothing" and "flattening" options for the Lambert Conformal projection not seen
    in the Cartopy version. By using Lambert Cylindrical over Lambert Conformal in Python,
    you will be able to create the "rectangular" style of coordinates not classically
    represented by a Lambert Conformal map. Additionally, Cartopy does not currently support
    adding tick marks to a projection like NCL, this is why these Python projections do not
    have this feature. The GeoCAT Team is actively adding to the list of convenience functions
    supported and hopes to add this functionality one day.
    '''

    # Plot data and create colorbar
    prec = t.plot.contourf(ax=ax,
                           cmap="Blues_r",
                           transform=ccrs.PlateCarree(),
                           levels=14,
                           add_colorbar=False)

    cbar_ticks = np.arange(0, 240, 20)
    cbar = plt.colorbar(prec,
                        orientation='horizontal',
                        pad=0.075,
                        shrink=0.8,
                        ticks=cbar_ticks)

    cbar.ax.tick_params(labelsize=10)
    plt.title(title, loc='center', y=1.17, size=15)
    plt.title(t.units, loc='right', y=1.08, size=14)
    plt.title("precipitation", loc='left', y=1.08, size=14)


Plot(
    2, 2, 1,
    ccrs.LambertConformal(central_longitude=45,
                          standard_parallels=(36, 55),
                          globe=ccrs.Globe()), "Lambert Conformal")
Plot(2, 2, 2, ccrs.LambertCylindrical(central_longitude=45),
     "Lambert Cylindrical")
Plot(2, 2, 3, ccrs.LambertAzimuthalEqualArea(central_longitude=45),
     "Lambert Azimuthal")