Plotting curvilinear with wrapping issue: suggestion on how to do it

[htonchia]

Description

Curvilinear grids are grids indexed by two dummy variables X and Y that are not defined in projection coordinates. The coordinates of the vertices can be given in latitudes longitudes but are not rectilinear in any projection. In some case, rectilinear data in a projection for which the parameters are missing and for which the vertices coordinates are known only in longitudes and latitudes can be handled as curvilinear grids. Currently Matplotlib and cartopy don't really handle the wrapping, or crossing of the projection boundary. There is (or will) be an exception. pcolormesh has a patch which (when updated #1418 and #1420) manages the overlapping cells produced by a wrapping grid.

The next figure shows a curvilinear dataset plotted with Matplolib. The first scatter plot shows the data and a third of the vertices. The second scatter plot shows the same data and vertices on a radial plot. The Atlantic Ocean is in the lower half center, with Groenland above, Europe on the right and North America on the right.

The contour, contourf, pcolor and pcolormesh plots are traversed by overlapping cells and contours.

When we plot the same data with cartopy in PlateCarree and Stereographic projections, we discover that the pcolormesh plots are fine while the pcolor, contour and contourf are traversed by overlapping data. The fine pcolormesh plots are due to the patch managing the overlapping data.

PLEASE NOTE, that this has been plotted with a geoaxes.py updated with #1418 and #1420 PRs. Without this update, the pcolormesh plots may not be so fine.

SUGGESTION I suggest to create a Geodetic_curvilinear CRS to handle those data. I suggest to create a function to handle overlaps in data in that Geodetic_curvilinear CRS, the same way it is handled in pcolormesh. When we get that Geodetic_curvilinear CRS in pcolor, contour or contourf, instead of rejecting it as it is not a projection, we would project the data in the axe.projection and handle the overlaps by masking the corresponding data vertices before sending the projected coordinates and the masked data to the Matplotib pcolor, contour or contourf function.

This should also fix issues #1400 #1151 (contourf case) #1225 #1076 and #522. SOLUTION BY PREPROCESSING THE DATA The overlaps cells can be masked before plotting the data with cartopy with pcolor, contour or contourf (even pcolormesh).

The result with the dataset above is the following:

The results with two of the datasets from issue #522 are:

Code to reproduce

import numpy as np
import numpy.ma as ma
import matplotlib.pyplot as plt
from netCDF4 import Dataset
import cartopy.crs as ccrs
try:
    import pykdtree.kdtree
    _is_pykdtree = True
except ImportError:
    import scipy.spatial
    _is_pykdtree = False

# file available here:
# https://github.com/htonchia/cartopy-plotting-examples/blob/master\
#/ice_conc_nh_polstere-100_multi_201912091200.nc

FILE = "ice_conc_nh_polstere-100_multi_201912091200.nc"
VARNAME = 'ice_conc'
XYNAME = 'xc', 'yc'
FIG = True

def data2curvi(xc, yc, zdata):
    """
    create a curvilinear dataset from data
    and returns the lat lon considering the target_xy in the proj:
        proj_data = ccrs.Stereographic(
        central_longitude=-45, central_latitude=90,
        true_scale_latitude=70)
    """
    nx = ny = 200

    Y = np.linspace(0, nx-1, nx)
    X = np.linspace(0, ny-1, ny)

    # from X Y curvi to stereopolaire target xc yc

    XX, YY = np.meshgrid(X, Y)
    theta = (YY + XX/4.) * np.pi / 100.
    radius = 10 * X + 0.1 * X**2

    # coordinates in km in graphic stereographic
    txc = radius * np.cos(theta)
    tyc = radius * np.sin(theta)

    target_xy = np.stack((txc.flatten(), tyc.flatten()), axis=-1)
    mxyc = np.meshgrid(xc, yc)
    xyc = np.stack((mxyc[0].flatten(), mxyc[1].flatten()), axis=-1)

    if _is_pykdtree:
        kdtree = pykdtree.kdtree.KDTree(xyc)
        dists, indices = kdtree.query(target_xy, k=1)
    else:
        # Versions of scipy >= v0.16 added the balanced_tree argument,
        # which caused the KDTree to hang with this input.
        try:
            kdtree = scipy.spatial.cKDTree(xyc, balanced_tree=False)
        except TypeError:
            kdtree = scipy.spatial.cKDTree(xyc)
        dists, indices = kdtree.query(target_xy, k=1)
    zdata1d = zdata.flatten()
    newdata = np.ma.array(zdata1d[indices])
    newdata = newdata.reshape(ny, nx)
    newdata = ma.masked_where(dists.reshape(ny,nx) > 100, newdata)

    proj_data = ccrs.Stereographic(
        central_longitude=-45, central_latitude=90,
        true_scale_latitude=70)

    lonlat = ccrs.Geodetic().transform_points(proj_data, txc * 1000, tyc * 1000)

    return lonlat[..., 0], lonlat[..., 1], newdata, XX, YY, txc, tyc

def main(file, varname, xy_name):
    proj_data = ccrs.Stereographic(
            central_longitude=-45, central_latitude=90,
            true_scale_latitude=70)

    dataset = Dataset(file)
    data = dataset.variables[varname][0,:,:]
    xc = dataset.variables[xy_name[0]][:]
    yc = dataset.variables[xy_name[1]][:]

    lons, lats, ZZ, X, Y, txc, tyc = data2curvi(xc, yc, zdata=data)

    plt.figure(figsize=(12, 6.4), dpi=100)
    plt.subplots_adjust(wspace=0.01, left=0, right=0.99)
    axe = plt.subplot(2, 3, 1)
    plt.scatter(lons, lats, 5, c=ZZ, cmap='bwr')
    plt.scatter(lons[::3,::3], lats[::3,::3], s=1, c=Y[::3,::3])
    axe.set_title('matplotlib scatter curvi')
    axe = plt.subplot(2, 3, 2)
    plt.pcolor(lons, lats, ZZ, cmap='bwr')
    axe.set_title('matplotlib pcolor curvi')
    axe = plt.subplot(2, 3, 3)
    plt.pcolormesh(lons, lats, ZZ, cmap='bwr')
    axe.set_title('matplotlib pcolormesh')
    axe = plt.subplot(2, 3, 4)
    plt.scatter(txc, tyc, c=ZZ, cmap='bwr')
    plt.scatter(txc[::3,::3], tyc[::3,::3], s=1, c=Y[::3,::3])
    axe.set_title('matplotlib radial plot curvi')

    axe = plt.subplot(2, 3, 5)
    plt.contour(lons, lats, ZZ, colors=['fuchsia'])
    axe.set_title('matplotlib contour curvi')
    axe = plt.subplot(2, 3, 6)
    plt.contourf(lons, lats, ZZ, cmap='bwr')
    axe.set_title('matplotlib contourf curvi')
    if FIG:
        forme = "img/curvi_matplotlib.png"
        plt.savefig(forme,
                        dpi=200, pad_inches=0)    

    range1=slice(None,None,None)
    range2=200
    #range1.extend(list(range(75,199)))

    transform = ccrs.PlateCarree(central_longitude=0.00001)

    plt.figure(figsize=(12, 6.4), dpi=100)
    plt.subplots_adjust(wspace=0.01, left=0, right=0.99)
    axe = plt.subplot(3, 3, 1, projection=ccrs.PlateCarree(central_longitude=0))
    axe.scatter(lons[range1,:range2], lats[range1,:range2], 1,c=ZZ[range1,:range2],
                transform=transform, cmap='bwr')
    axe.set_title('cartopy scatter curvi PC')
    axe.coastlines()
    axe = plt.subplot(3, 3, 2, projection=ccrs.PlateCarree(central_longitude=0))
    axe.pcolor(lons[range1,:range2], lats[range1,:range2], ZZ[range1,:range2],
               cmap='bwr', transform=transform)
    axe.set_title('cartopy pcolor curvi PC')
    axe.coastlines()
    axe = plt.subplot(3, 3, 3, projection=ccrs.PlateCarree(central_longitude=0))
    axe.pcolormesh(lons[range1,:range2], lats[range1,:range2], ZZ[range1,:range2],
                   cmap='bwr', transform=transform)
    axe.set_title('cartopy pcolormesh curvi PC')
    axe.coastlines()
    axe = plt.subplot(3, 3, 4, projection=ccrs.PlateCarree(central_longitude=0))
    axe.contour(lons[range1,:range2], lats[range1,:range2], ZZ[range1,:range2],
                colors=['fuchsia'], transform=transform)
    axe.set_title('cartopy contour curvi PC')
    axe.coastlines()
    axe = plt.subplot(3, 3, 5, projection=ccrs.PlateCarree(central_longitude=0))
    axe.contourf(lons[range1,:range2], lats[range1,:range2], ZZ[range1,:range2],
                 cmap='bwr', transform=transform)
    axe.set_title('cartopy contourf curvi PC')
    axe.coastlines()
    axe = plt.subplot(3, 3, 6, projection=proj_data)
    axe.pcolor(lons[range1,:range2], lats[range1,:range2], ZZ[range1,:range2],
               cmap='bwr', transform=transform)
    axe.set_title('cartopy pcolor curvi PC->Stere')
    axe.coastlines()
    axe = plt.subplot(3, 3, 7, projection=proj_data)
    axe.pcolormesh(lons[range1,:range2], lats[range1,:range2], ZZ[range1,:range2],
                   cmap='bwr', transform=transform)
    axe.set_title('cartopy pcolormesh curvi PC->Stere')
    axe.coastlines()
    axe = plt.subplot(3, 3, 8, projection=proj_data)
    axe.contour(lons[range1,:range2], lats[range1,:range2], ZZ[range1,:range2],
                colors=['fuchsia'], transform=transform)
    axe.set_title('cartopy contour curvi PC->Stere')
    axe.coastlines()
    axe = plt.subplot(3, 3, 9, projection=proj_data)
    axe.contourf(lons[range1,:range2], lats[range1,:range2], ZZ[range1,:range2],
                 cmap='bwr', transform=transform)
    axe.set_title('cartopy contourf curvi PC->Stere')
    axe.coastlines()
    if FIG:
        forme = "img/curvi_cartopy.png"
        plt.savefig(forme,
                        dpi=200, pad_inches=0)

if __name__ == '__main__':
    main(file=FILE, varname=VARNAME, xy_name=XYNAME)

CODE WITH THE PREPROCESSING FIX FUNCTION z_masked_overlap

import numpy as np
import numpy.ma as ma
import matplotlib.pyplot as plt
from netCDF4 import Dataset
import cartopy.crs as ccrs
try:
    import pykdtree.kdtree
    _IS_PYKDTREE = True
except ImportError:
    import scipy.spatial
    _IS_PYKDTREE = False

# file available here:
# https://github.com/htonchia/cartopy-plotting-examples/blob/master\
#/ice_conc_nh_polstere-100_multi_201912091200.nc

FILE = "ice_conc_nh_polstere-100_multi_201912091200.nc"
VARNAME = 'ice_conc'
XYNAME = 'xc', 'yc'

FILE2 = "https://vesg.ipsl.upmc.fr/thredds/dodsC/IPSLFS/brocksce/matplotlib/ORCA2.0_ex1.nc"
FILE2 = "tos.nc"
VARNAME2 = 'tos'
XYNAME2 = 'lon', 'lat'

FILE3 = "https://vesg.ipsl.upmc.fr/thredds/dodsC/IPSLFS/brocksce/matplotlib/ORCA0.25_ex1.nc"
FILE3 = "var.nc"
VARNAME3 = 'VAR'
XYNAME3 = 'NAV_LON', 'NAV_LAT'

FIG = True

def data2curvi(xc, yc, zdata, proj_data):
    """
    create a curvilinear dataset from data
    and gives the lat lon considering the target_xy in the proj:
        proj_data = ccrs.Stereographic(
        central_longitude=-45, central_latitude=90,
        true_scale_latitude=70)
    """
    nx = ny = 200

    Y = np.linspace(0, ny-1, ny)
    X = np.linspace(0, nx-1, nx)

    # from X Y curvi to stereopolaire target xc yc

    XX, YY = np.meshgrid(X, Y)
    theta = (YY + XX/4.) * np.pi / 100.
    radius = 10 * X + 0.1 * X**2

    # coordinates in km in graphic stereographic
    txc = radius * np.cos(theta)
    tyc = radius * np.sin(theta)

    target_xy = np.stack((txc.flatten(), tyc.flatten()), axis=-1)
    mxyc = np.meshgrid(xc, yc)
    xyc = np.stack((mxyc[0].flatten(), mxyc[1].flatten()), axis=-1)

    if _IS_PYKDTREE:
        kdtree = pykdtree.kdtree.KDTree(xyc)
        dists, indices = kdtree.query(target_xy, k=1)
    else:
        # Versions of scipy >= v0.16 added the balanced_tree argument,
        # which caused the KDTree to hang with this input.
        try:
            kdtree = scipy.spatial.cKDTree(xyc, balanced_tree=False)
        except TypeError:
            kdtree = scipy.spatial.cKDTree(xyc)
        dists, indices = kdtree.query(target_xy, k=1)
    zdata1d = zdata.flatten()
    newdata = np.ma.array(zdata1d[indices])
    newdata = newdata.reshape(ny, nx)
    newdata = ma.masked_where(dists.reshape(ny, nx) > 100, newdata)

    lonlat = ccrs.Geodetic().transform_points(proj_data, txc * 1000, tyc * 1000)

    return lonlat[..., 0], lonlat[..., 1], newdata, XX, YY, txc, tyc

def z_masked_overlap(axe, X, Y, Z, source_projection=None):
    """
    for data in projection axe.projection
    find and mask the overlaps (more 1/2 the axe.projection range)

    X, Y either the coordinates in axe.projection or longitudes latitudes
    Z the data
    operation one of 'pcorlor', 'pcolormesh', 'countour', 'countourf'

    if source_projection is a geodetic CRS data is in geodetic coordinates
    and should first be projected in axe.projection

    X, Y are 2D same dimension as Z for contour and contourf
    same dimension as Z or with an extra row and column for pcolor
    and pcolormesh

    return ptx, pty, Z
    """
    if not hasattr(axe, 'projection'):
        return X, Y, Z
    if not isinstance(axe.projection, ccrs.Projection):
        return X, Y, Z

    if len(X.shape) != 2 or len(Y.shape) != 2:
        return X, Y, Z

    if (source_projection is not None and
            isinstance(source_projection, ccrs.Geodetic)):
        transformed_pts = axe.projection.transform_points(
            source_projection, X, Y)
        ptx, pty = transformed_pts[..., 0], transformed_pts[..., 1]
    else:
        ptx, pty = X, Y

    with np.errstate(invalid='ignore'):
        # diagonals have one less row and one less columns
        diagonal0_lengths = np.hypot(
            ptx[1:, 1:] - ptx[:-1, :-1],
            pty[1:, 1:] - pty[:-1, :-1]
        )
        diagonal1_lengths = np.hypot(
            ptx[1:, :-1] - ptx[:-1, 1:],
            pty[1:, :-1] - pty[:-1, 1:]
        )
        to_mask = (
            (diagonal0_lengths > (
                abs(axe.projection.x_limits[1]
                    - axe.projection.x_limits[0])) / 2) |
            np.isnan(diagonal0_lengths) |
            (diagonal1_lengths > (
                abs(axe.projection.x_limits[1]
                    - axe.projection.x_limits[0])) / 2) |
            np.isnan(diagonal1_lengths)
        )

        # TODO check if we need to do something about surrounding vertices

        # add one extra colum and row for contour and contourf
        if (to_mask.shape[0] == Z.shape[0] - 1 and
                to_mask.shape[1] == Z.shape[1] - 1):
            to_mask_extended = np.zeros(Z.shape, dtype=bool)
            to_mask_extended[:-1, :-1] = to_mask
            to_mask_extended[-1, :] = to_mask_extended[-2, :]
            to_mask_extended[:, -1] = to_mask_extended[:, -2]
            to_mask = to_mask_extended
        if np.any(to_mask):

            Z_mask = getattr(Z, 'mask', None)
            to_mask = to_mask if Z_mask is None else to_mask | Z_mask

            Z = ma.masked_where(to_mask, Z)

        return ptx, pty, Z

def main(file, varname, xy_name, to_curvi_model=True,
         plotkwargs=None):
    """
    file: nc file
    varname: the variable to plot
    xy_name: the longitude and latitudes variables
    to_curvi_model: transform data to a curvilenar dataset
    if false
    assume lons and lats are 2D
    """
    dataset = Dataset(file)
    data = dataset.variables[varname][0, :, :]
    xc = dataset.variables[xy_name[0]][:]
    yc = dataset.variables[xy_name[1]][:]
    proj_data = ccrs.Stereographic(
        central_longitude=-45, central_latitude=90,
        true_scale_latitude=70)

    if to_curvi_model:
        lons, lats, ZZ, X, Y, txc, tyc = data2curvi(
            xc, yc, zdata=data, proj_data=proj_data)
    else:
        lons, lats, ZZ = xc, yc, data
        nnj, nni = lons.shape
        X, Y = np.meshgrid(
            np.linspace(0, nni - 1, nni),
            np.linspace(0, nnj - 1, nnj))

    operations = ['pcolor', 'pcolormesh', 'contour', 'contourf']

    plotkwargs = {
        'pcolor' : {'cmap' : 'bwr'},
        'pcolormesh' : {'cmap' : 'bwr'},
        'contour' : {'colors' : ['fuchsia']},
        'contourf' : {'cmap' : 'bwr'}
        } if plotkwargs is None else plotkwargs

    projections = [
        ccrs.PlateCarree(central_longitude=0),
        proj_data,
        ccrs.Robinson(central_longitude=180)]

    if to_curvi_model:
        cols = 3
        rows = 2
        plt.figure(figsize=(12, 6.4), dpi=100)
        plt.subplots_adjust(wspace=0.01, left=0, right=0.99)
        plt.suptitle('{} Matplotlib only'.format(varname))

        # showing the curvilinear grid and data
        axe = plt.subplot(rows, cols, 1)
        plt.scatter(lons, lats, 5, c=ZZ, cmap=plotkwargs['pcolor']['cmap'])
        plt.scatter(lons[::3, ::3], lats[::3, ::3], s=1, c=Y[::3, ::3])
        axe.set_title('matplotlib scatter curvi')

        axe = plt.subplot(rows, cols, 4)
        plt.scatter(txc, tyc, c=ZZ, cmap=plotkwargs['pcolor']['cmap'])
        plt.scatter(txc[::3, ::3], tyc[::3, ::3], s=1, c=Y[::3, ::3])
        axe.set_title('matplotlib radial plot curvi')

        # maplotlib plots

        for iop, operation in enumerate(operations):
            iii = iop + 2 if iop < 2 else iop + 3
            axe = plt.subplot(2, 3, iii)
            getattr(axe, operation)(lons, lats, ZZ,
                                    **plotkwargs[operation])
            axe.set_title('matplotlib {} curvi'.format(operation))

        if FIG:
            forme = "img/curvi_{}_matplotlib.png"
            plt.savefig(forme.format(varname),
                        dpi=200, pad_inches=0)

    #range1.extend(list(range(75,199)))

    transform = ccrs.PlateCarree(central_longitude=0.00001)
    # a trick to force going into pcolormesh patch

    plt.figure(figsize=(12, 6.4), dpi=100)
    plt.subplots_adjust(wspace=0.01, left=0, right=0.99)
    plt.suptitle('{} Cartopy with pcolormesh bug fix'.format(varname))
    cols = len(operations)
    rows = len(projections) + 1
    #scatter
    axe = plt.subplot(rows, cols, 1, projection=ccrs.PlateCarree(central_longitude=0))
    kwargs = {'cmap' : plotkwargs['pcolor']['cmap']}
    if 'vmin' in plotkwargs['pcolor'] and 'vmax' in plotkwargs['pcolor']:
        kwargs['vmin'] = plotkwargs['pcolor']['vmin']
        kwargs['vmax'] = plotkwargs['pcolor']['vmax']
    axe.scatter(lons, lats, 1, c=ZZ,
                transform=transform, **kwargs)
    axe.set_title('cartopy scatter curvi PC')
    axe.coastlines()

    iii = cols
    for projection in projections:
        for operation in operations:
            iii += 1
            axe = plt.subplot(rows, cols, iii, projection=projection)
            getattr(axe, operation)(lons, lats, ZZ,
                                    **plotkwargs[operation],
                                    transform=transform)
            axe.set_title('cartopy {} curvi {}'.format(
                operation,
                projection.proj4_params['proj']))
            axe.coastlines()

    if FIG:
        forme = "img/curvi_{}_cartopy.png"
        plt.savefig(forme.format(varname),
                    dpi=200, pad_inches=0)

    plt.figure(figsize=(12, 6.4), dpi=100)
    plt.subplots_adjust(wspace=0.01, left=0, right=0.99)
    plt.suptitle('{} Cartopy with pcolormesh bug fix and overlap preprocess'.
                 format(varname))
    rows = rows - 1
#    axe = plt.subplot(rows, cols, 1, projection=ccrs.PlateCarree(central_longitude=0))
#    axe.scatter(lons, lats, 1, c=ZZ,
#                transform=transform, cmap='bwr')
#    axe.set_title('cartopy scatter curvi PC')
#    axe.coastlines()

    iii = 0

    for projection in projections:
        for operation in operations:
            iii += 1
            axe = plt.subplot(rows, cols, iii, projection=projection)
            X, Y, maskedZ = z_masked_overlap(
                axe, lons, lats, ZZ,
                source_projection=ccrs.Geodetic())
            getattr(axe, operation)(X, Y, maskedZ,
                                    **plotkwargs[operation])
            axe.set_title('cartopy {} curvi {}'.format(
                operation,
                projection.proj4_params['proj']))
            axe.coastlines()

    if FIG:
        forme = "img/curvi_{}_cartopy_FIX.png"
        plt.savefig(forme.format(varname),
                    dpi=200, pad_inches=0)

if __name__ == '__main__':
    cmap = 'bwr'
    colors = ['fuchsia']
    pkwargs = {
        'pcolor' : {'cmap': cmap},
        'pcolormesh' : {'cmap': cmap},
        'contour' : {'colors': colors},
        'contourf' : {'cmap': cmap}
        }
    main(
        file=FILE,
        varname=VARNAME,
        xy_name=XYNAME,
        to_curvi_model=True,
        plotkwargs=pkwargs
        )
    cmap = 'plasma'
    pkwargs = {
        'pcolor' : {'cmap' : cmap},
        'pcolormesh' : {'cmap' : cmap},
        'contour' : {'colors' : colors},
        'contourf' : {'cmap' : cmap}
        }

    main(
        file=FILE2,
        varname=VARNAME2,
        xy_name=XYNAME2,
        to_curvi_model=False,
        plotkwargs=pkwargs)
    levels = np.arange(0, 60e-9, 5e-9)
    pkwargs = {
        'pcolor' : {'cmap' : cmap, 'vmin' : 0, 'vmax' : 60e-9},
        'pcolormesh' : {'cmap' : cmap, 'vmin' : 0, 'vmax' : 60e-9},
        'contour' : {'colors' : colors, 'levels' : levels},
        'contourf' : {'cmap' : cmap, 'levels' : levels}
        }
    # BEWARE this part is quite time consuming
    main(
        file=FILE3,
        varname=VARNAME3,
        xy_name=XYNAME3,
        to_curvi_model=False,
        plotkwargs=pkwargs)

Traceback

Full environment definition

### Operating system ### Cartopy version ### conda list ``` ``` ### pip list ``` ```

[jonnyhtw]

Hi there,

Is there a workaround for this as yet? I'm trying to plot contours on top of a `pcolormesh' instance with cartopy but I consistently get the wrapping issue. The ocean field I am trying to plot is regional, not global if that makes any difference?

Thanks!