Curvilinear Datagrid is not displayed correctly

[davidmnielsen]

The datagrid is not displayed correctly, in this specific case for the curvilinear grid of MPIOM when trying the datagrid option at psy.plot.mapplot:

maps = psy.plot.mapplot(file_lr,
                        name='area',
                        cmap='Blues',
                        maskleq=0,
                        stock_img=True,
                        lonlatbox=[-180, 180, 60, 90],
                        extend='max',
                        load=True,
                        projection='northpole',
                        xgrid=False, ygrid=False,
                        clabel='{desc}',title='LR',
                        bounds=np.arange(0,10000000000,1000000000),
                        datagrid=dict(color='k', linewidth=0.1),
                        tight=True)

A full reproducible example is available at Github, with all the auxiliary files needed: https://github.com/davidmnielsen/gridcell_area

[Chilipp]

Hi @davidmnielsen! Thanks for reporting this :smile: As already mentioned per mail, a workaround is to use matplotlibs edgecolor property via

maps = psy.plot.mapplot(file_lr,
                        name='area',
                        cmap='Blues',
                        maskleq=0,
                        stock_img=True,
                        lonlatbox=[-180, 180, 60, 90],
                        extend='max',
                        load=True,
                        projection='northpole',
                        xgrid=False, ygrid=False,
                        clabel='{desc}',title='LR',
                        bounds=np.arange(0,10000000000,1000000000),
                        datagrid=dict(color='k', linewidth=0.1),
                        tight=True,
                        post="self.plotter.plot.mappable.set_edgecolor('k')")

But I am working on it. The problem appears only for grid with 2D coordinates, such as the curvilinear grid you are using here.

[Chilipp]

Hey! It seems like I have to fix an old issue with the visualization of curvilinear grids in order to solve this problem.

The point is: Matplotlib actually visualizes your data wrong. When typing

plot.pcolormesh(x, y, data)

with data.shape == y.shape (as it is the case for curvilinear grids), matplotlib gets rid of the last row and last column of the data and handles the grid cell centers in y as the cell node coordinates (see matplotlibs source code). However, as I mentioned already in the circumpolar grid example, this introduces an error in the visualization.

psyplot usually interpolates the bounds to avoid this error, but I did not yet find a good solution on how to do this with 2D coordinates. So far, I disable the interpolation with the interp_bounds formatoption, and, if desired, one can interpolate the bounds using scipys interp2d function (see psyplot.data.CFDecoder._infer_interval_breaks). This function however does not work with NaN in the data and as such, it fails at the boundaries of the curvilinear grid.

Long story short: If there is no suggestion on how to get the grid cell boundaries out of the 2D coordinates, I have to disable the visualization of data grids for curvilinear data. Any idea how to solve this?

[davidmnielsen]

Hi!

I found a solution, or at least for the visualization outcome I was aiming for. Somehow (extensive Stack Exchange digging), I got pcolormesh to work in Cartopy.

Thank you a lot for your time and effort to help us! :)

def circlemap():
    import numpy as np
    import matplotlib.path as mpath
    theta = np.linspace(0, 2*np.pi, 100)
    center, radius = [0.5, 0.5], 0.5
    verts = np.vstack([np.sin(theta), np.cos(theta)]).T
    circle = mpath.Path(verts * radius + center)
    return circle

fig = plt.figure(figsize=(8, 8))
ax = fig.add_subplot(111,projection=ccrs.NorthPolarStereo())
ax.set_extent([-180, 180, 63, 90], crs=ccrs.PlateCarree())
ax.coastlines()
circle=circlemap()
ax.set_boundary(circle, transform=ax.transAxes)

mesh = ax.pcolormesh(lon, lat, area*factor, cmap='Blues',
                     transform=ccrs.PlateCarree(),
                     vmin=0, vmax=10000,
                     antialiased=True, edgecolor='k', linewidth=0.3)
plt.colorbar(mesh, orientation='horizontal', shrink=0.75,label='grid-cell area [km^2]')

plt.show()

I added a full reproducible Jupyter-Notebook to the repository: https://github.com/davidmnielsen/gridcell_area

[davidmnielsen]

Oh, and yes, adding the argument post="self.plotter.plot.mappable.set_edgecolor('k')" indeed does the trick in psyplot.

[Chilipp]

Alright! Thanks for the confirmation :) I'll decide on the next few days how to handle the Datagrid for Curvilinear grids and then close this issue. Any input or suggestions are welcomed!