SciTools / cartopy

Cartopy - a cartographic python library with matplotlib support
https://scitools.org.uk/cartopy/docs/latest
BSD 3-Clause "New" or "Revised" License
1.41k stars 359 forks source link

new shading default can mess up plots #1638

Closed mathause closed 2 years ago

mathause commented 4 years ago

Description

The new shading default https://github.com/SciTools/cartopy/blob/08db54f2acbf216d19c3d8245dd375fc9dcf77ba/lib/cartopy/mpl/geoaxes.py#L1669 can lead to artifacts in the plot which do not happen with shading="flat". This happens for coordinates that cross the zero line. As far as I can see matplotlib plans to entirely remove this option?

Code to reproduce

I was too lazy to create an example without xarray as dependency - let me know if you need one:

Artifacts:

import cartopy.crs as ccrs
import xarray as xr
rasm = xr.tutorial.load_dataset("rasm")

# choose a projection
proj = ccrs.NorthPolarStereo()

ax = plt.subplot(111, projection=proj)
ax.set_global()

rasm.isel(time=1).Tair.plot.pcolormesh(
    ax=ax, x="xc", y="yc", transform=ccrs.PlateCarree()
)

ax.set_extent([-180, 180, 43, 90], ccrs.PlateCarree())

ax.coastlines();

shading_auto

No artifacts:

f = plt.fugure()

ax = plt.subplot(111, projection=proj)
ax.set_global()

rasm.isel(time=1).Tair.plot.pcolormesh(
    ax=ax, x="xc", y="yc", transform=ccrs.PlateCarree(), shading="flat"
)

ax.set_extent([-180, 180, 43, 90], ccrs.PlateCarree())

ax.coastlines();

shading_flat

Full environment definition ### Operating system linux ### Cartopy version 0.18.0 ### conda list ``` # packages in environment at /home/mathause/conda/envs/xarray_dev: # # Name Version Build Channel _libgcc_mutex 0.1 conda_forge conda-forge _openmp_mutex 4.5 1_gnu conda-forge affine 2.3.0 py_0 conda-forge antlr-python-runtime 4.7.2 py38_1001 conda-forge apipkg 1.5 py_0 conda-forge appdirs 1.4.3 py_1 conda-forge asciitree 0.3.3 py_2 conda-forge attrs 20.1.0 pyh9f0ad1d_0 conda-forge backcall 0.2.0 pyh9f0ad1d_0 conda-forge backports 1.0 py_2 conda-forge backports.functools_lru_cache 1.6.1 py_0 conda-forge beautifulsoup4 4.9.1 py_1 conda-forge black 19.10b0 py_4 conda-forge bokeh 2.1.1 py38h32f6830_0 conda-forge boost-cpp 1.72.0 h8e57a91_0 conda-forge boto3 1.14.47 pyh9f0ad1d_0 conda-forge botocore 1.17.47 pyh9f0ad1d_0 conda-forge bottleneck 1.3.2 py38h8790de6_1 conda-forge brotlipy 0.7.0 py38h1e0a361_1000 conda-forge bzip2 1.0.8 h516909a_2 conda-forge c-ares 1.16.1 h516909a_0 conda-forge ca-certificates 2020.6.20 hecda079_0 conda-forge cairo 1.16.0 hcf35c78_1003 conda-forge cartopy 0.18.0 py38h172510d_0 conda-forge cdat_info 8.2.1 pyh9f0ad1d_1 conda-forge cdms2 3.1.5 pypi_0 pypi cdtime 3.1.4 py38hac60b08_0 conda-forge certifi 2020.6.20 py38h32f6830_0 conda-forge cf-units 2.1.4 py38h8790de6_0 conda-forge cffi 1.14.1 py38h5bae8af_0 conda-forge cfgrib 0.9.8.4 py_0 conda-forge cfitsio 3.470 hce51eda_6 conda-forge cftime 1.2.1 py38h8790de6_0 conda-forge chardet 3.0.4 py38h32f6830_1006 conda-forge click 7.1.2 pyh9f0ad1d_0 conda-forge click-plugins 1.1.1 py_0 conda-forge cligj 0.5.0 py_0 conda-forge cloudpickle 1.5.0 py_0 conda-forge coverage 5.2.1 py38h1e0a361_0 conda-forge coveralls 2.1.2 py_0 conda-forge cryptography 3.0 py38h766eaa4_0 conda-forge curl 7.71.1 he644dc0_5 conda-forge cycler 0.10.0 py_2 conda-forge cytoolz 0.10.1 py38h516909a_0 conda-forge dask 2.23.0 py_0 conda-forge dask-core 2.23.0 py_0 conda-forge decorator 4.4.2 py_0 conda-forge distarray 2.12.2 py_1 conda-forge distributed 2.23.0 py38h32f6830_0 conda-forge docopt 0.6.2 py_1 conda-forge docutils 0.15.2 py38_0 conda-forge eccodes 2.17.0 h59f7be3_1 conda-forge esmf 8.0.0 nompi_hb0fcdcb_6 conda-forge esmpy 8.0.0 nompi_py38hf0e99fa_1 conda-forge execnet 1.7.1 py_0 conda-forge expat 2.2.9 he1b5a44_2 conda-forge fasteners 0.14.1 py_3 conda-forge flake8 3.8.3 py_1 conda-forge fontconfig 2.13.1 h86ecdb6_1001 conda-forge freetype 2.10.2 he06d7ca_0 conda-forge freexl 1.0.5 h516909a_1002 conda-forge fsspec 0.8.0 py_0 conda-forge future 0.18.2 py38h32f6830_1 conda-forge g2clib 1.6.0 hf3f1b0b_9 conda-forge gdal 3.0.4 py38h172510d_6 conda-forge geos 3.8.1 he1b5a44_0 conda-forge geotiff 1.5.1 h05acad5_10 conda-forge gettext 0.19.8.1 hc5be6a0_1002 conda-forge giflib 5.2.1 h516909a_2 conda-forge glib 2.65.0 h6f030ca_0 conda-forge h5netcdf 0.8.1 py_0 conda-forge h5py 2.10.0 nompi_py38h513d04c_102 conda-forge hdf4 4.2.13 hf30be14_1003 conda-forge hdf5 1.10.5 nompi_h3c11f04_1104 conda-forge hdfeos2 2.20 h64bfcee_1000 conda-forge hdfeos5 5.1.16 h8b6279f_5 conda-forge heapdict 1.0.1 py_0 conda-forge hypothesis 5.27.0 py_0 conda-forge icu 64.2 he1b5a44_1 conda-forge idna 2.10 pyh9f0ad1d_0 conda-forge importlib-metadata 1.7.0 py38h32f6830_0 conda-forge importlib_metadata 1.7.0 0 conda-forge iniconfig 1.0.1 pyh9f0ad1d_0 conda-forge ipython 7.17.0 py38h1cdfbd6_0 conda-forge ipython_genutils 0.2.0 py_1 conda-forge iris 2.4.0 py38_0 conda-forge isort 5.4.2 py38h32f6830_0 conda-forge jasper 1.900.1 h07fcdf6_1006 conda-forge jedi 0.17.2 py38h32f6830_0 conda-forge jinja2 2.11.2 pyh9f0ad1d_0 conda-forge jmespath 0.10.0 pyh9f0ad1d_0 conda-forge jpeg 9d h516909a_0 conda-forge json-c 0.13.1 hbfbb72e_1002 conda-forge jsonschema 3.2.0 py38h32f6830_1 conda-forge jupyter_core 4.6.3 py38h32f6830_1 conda-forge kealib 1.4.13 hec59c27_0 conda-forge kiwisolver 1.2.0 py38hbf85e49_0 conda-forge krb5 1.17.1 hfafb76e_2 conda-forge lazy-object-proxy 1.5.1 py38h1e0a361_0 conda-forge lcms2 2.11 hbd6801e_0 conda-forge ld_impl_linux-64 2.34 hc38a660_9 conda-forge libaec 1.0.4 he1b5a44_1 conda-forge libblas 3.8.0 17_openblas conda-forge libcblas 3.8.0 17_openblas conda-forge libcdms 3.1.2 hf60d256_111 conda-forge libcf 1.0.3 py38h2f41aa0_108 conda-forge libcurl 7.71.1 hcdd3856_5 conda-forge libdap4 3.20.6 h1d1bd15_1 conda-forge libdrs 3.1.2 hc2e2db3_112 conda-forge libdrs_f 3.1.2 hae7e664_110 conda-forge libedit 3.1.20191231 he28a2e2_2 conda-forge libev 4.33 h516909a_0 conda-forge libffi 3.2.1 he1b5a44_1007 conda-forge libgcc-ng 9.3.0 h24d8f2e_15 conda-forge libgdal 3.0.4 h3dfc09a_6 conda-forge libgfortran-ng 7.5.0 hdf63c60_15 conda-forge libgomp 9.3.0 h24d8f2e_15 conda-forge libiconv 1.16 h516909a_0 conda-forge libkml 1.3.0 hb574062_1011 conda-forge liblapack 3.8.0 17_openblas conda-forge libllvm10 10.0.1 he513fc3_1 conda-forge libnetcdf 4.7.4 nompi_h9f9fd6a_101 conda-forge libnghttp2 1.41.0 hab1572f_1 conda-forge libopenblas 0.3.10 pthreads_hb3c22a3_4 conda-forge libpng 1.6.37 hed695b0_2 conda-forge libpq 12.3 h5513abc_0 conda-forge libspatialite 4.3.0a h2482549_1038 conda-forge libssh2 1.9.0 hab1572f_5 conda-forge libstdcxx-ng 9.3.0 hdf63c60_15 conda-forge libtiff 4.1.0 hc7e4089_6 conda-forge libuuid 2.32.1 h14c3975_1000 conda-forge libwebp-base 1.1.0 h516909a_3 conda-forge libxcb 1.13 h14c3975_1002 conda-forge libxml2 2.9.10 hee79883_0 conda-forge libxslt 1.1.33 h31b3aaa_0 conda-forge llvmlite 0.34.0 py38h4f45e52_0 conda-forge locket 0.2.0 py_2 conda-forge lxml 4.5.2 py38hbb43d70_0 conda-forge lz4-c 1.9.2 he1b5a44_3 conda-forge markupsafe 1.1.1 py38h1e0a361_1 conda-forge matplotlib 3.3.1 1 conda-forge matplotlib-base 3.3.1 py38h91b0d89_1 conda-forge mccabe 0.6.1 py_1 conda-forge mechanicalsoup 0.12.0 py_0 conda-forge monotonic 1.5 py_0 conda-forge more-itertools 8.4.0 py_0 conda-forge msgpack-python 1.0.0 py38hbf85e49_1 conda-forge mypy 0.780 py_0 conda-forge mypy_extensions 0.4.3 py38h32f6830_1 conda-forge nbformat 5.0.7 py_0 conda-forge nc-time-axis 1.2.0 py_1 conda-forge ncurses 6.2 he1b5a44_1 conda-forge netcdf-fortran 4.5.2 nompi_h45d7149_104 conda-forge netcdf4 1.5.3 nompi_py38heb6102f_103 conda-forge numba 0.51.0 py38hc5bc63f_0 conda-forge numbagg 0.1 pypi_0 pypi numcodecs 0.6.4 py38he1b5a44_0 conda-forge numpy 1.19.1 py38hbc27379_2 conda-forge olefile 0.46 py_0 conda-forge openblas 0.3.10 pthreads_hf183345_4 conda-forge openjpeg 2.3.1 h981e76c_3 conda-forge openssl 1.1.1g h516909a_1 conda-forge owslib 0.20.0 py_0 conda-forge packaging 20.4 pyh9f0ad1d_0 conda-forge pandas 1.1.0 py38h950e882_0 conda-forge parso 0.7.1 pyh9f0ad1d_0 conda-forge partd 1.1.0 py_0 conda-forge pathspec 0.8.0 pyh9f0ad1d_0 conda-forge patsy 0.5.1 py_0 conda-forge pcre 8.44 he1b5a44_0 conda-forge pexpect 4.8.0 py38h32f6830_1 conda-forge pickleshare 0.7.5 py38h32f6830_1001 conda-forge pillow 7.2.0 py38h9776b28_1 conda-forge pint 0.14 py_0 conda-forge pip 20.2.2 py_0 conda-forge pixman 0.38.0 h516909a_1003 conda-forge pluggy 0.13.1 py38h32f6830_2 conda-forge poppler 0.67.0 h14e79db_8 conda-forge poppler-data 0.4.9 1 conda-forge postgresql 12.3 h8573dbc_0 conda-forge proj 7.0.0 h966b41f_5 conda-forge prompt-toolkit 3.0.6 py_0 conda-forge pseudonetcdf 3.1.0 py_1 conda-forge psutil 5.7.2 py38h1e0a361_0 conda-forge pthread-stubs 0.4 h14c3975_1001 conda-forge ptyprocess 0.6.0 py_1001 conda-forge py 1.9.0 pyh9f0ad1d_0 conda-forge pycodestyle 2.6.0 pyh9f0ad1d_0 conda-forge pycparser 2.20 pyh9f0ad1d_2 conda-forge pydap 3.2.2 py38_1000 conda-forge pyepsg 0.4.0 py_0 conda-forge pyflakes 2.2.0 pyh9f0ad1d_0 conda-forge pygments 2.6.1 py_0 conda-forge pyke 1.1.1 py38h32f6830_1002 conda-forge pynio 1.5.5 py38h031d99c_12 conda-forge pyopenssl 19.1.0 py_1 conda-forge pyparsing 2.4.7 pyh9f0ad1d_0 conda-forge pyproj 2.6.1.post1 py38h7521cb9_0 conda-forge pyrsistent 0.16.0 py38h1e0a361_0 conda-forge pyshp 2.1.0 py_0 conda-forge pysocks 1.7.1 py38h32f6830_1 conda-forge pytest 6.0.1 py38h32f6830_0 conda-forge pytest-cov 2.10.1 pyh9f0ad1d_0 conda-forge pytest-env 0.6.2 py_0 conda-forge pytest-forked 1.2.0 pyh9f0ad1d_0 conda-forge pytest-xdist 2.0.0 py_0 conda-forge python 3.8.5 h1103e12_4_cpython conda-forge python-dateutil 2.8.1 py_0 conda-forge python_abi 3.8 1_cp38 conda-forge pytz 2020.1 pyh9f0ad1d_0 conda-forge pyyaml 5.3.1 py38h1e0a361_0 conda-forge rasterio 1.1.5 py38h033e0f6_1 conda-forge readline 8.0 he28a2e2_2 conda-forge regex 2020.7.14 py38h1e0a361_0 conda-forge regrid2 3.1.5 pypi_0 pypi requests 2.24.0 pyh9f0ad1d_0 conda-forge s3transfer 0.3.3 py38h32f6830_1 conda-forge scipy 1.5.2 py38h8c5af15_0 conda-forge seaborn 0.10.1 1 conda-forge seaborn-base 0.10.1 py_1 conda-forge setuptools 49.6.0 py38h32f6830_0 conda-forge shapely 1.7.1 py38hc7361b7_0 conda-forge six 1.15.0 pyh9f0ad1d_0 conda-forge snuggs 1.4.7 py_0 conda-forge sortedcontainers 2.2.2 pyh9f0ad1d_0 conda-forge soupsieve 2.0.1 py_1 conda-forge sparse 0.11.0 py_0 conda-forge sqlite 3.33.0 h4cf870e_0 conda-forge statsmodels 0.11.1 py38h1e0a361_2 conda-forge tbb 2020.1 hc9558a2_0 conda-forge tblib 1.6.0 py_0 conda-forge tiledb 1.7.7 h8efa9f0_3 conda-forge tk 8.6.10 hed695b0_0 conda-forge toml 0.10.1 pyh9f0ad1d_0 conda-forge toolz 0.10.0 py_0 conda-forge tornado 6.0.4 py38h1e0a361_1 conda-forge traitlets 4.3.3 py38h32f6830_1 conda-forge typed-ast 1.4.1 py38h516909a_0 conda-forge typing_extensions 3.7.4.2 py_0 conda-forge tzcode 2020a h516909a_0 conda-forge udunits2 2.2.27.6 h4e0c4b3_1001 conda-forge urllib3 1.25.10 py_0 conda-forge wcwidth 0.2.5 pyh9f0ad1d_1 conda-forge webob 1.8.6 py_0 conda-forge wheel 0.35.1 pyh9f0ad1d_0 conda-forge xarray 0.16.1.dev53+g43a2a4bd dev_0 xerces-c 3.2.2 h8412b87_1004 conda-forge xorg-kbproto 1.0.7 h14c3975_1002 conda-forge xorg-libice 1.0.10 h516909a_0 conda-forge xorg-libsm 1.2.3 h84519dc_1000 conda-forge xorg-libx11 1.6.11 h516909a_0 conda-forge xorg-libxau 1.0.9 h14c3975_0 conda-forge xorg-libxdmcp 1.1.3 h516909a_0 conda-forge xorg-libxext 1.3.4 h516909a_0 conda-forge xorg-libxrender 0.9.10 h516909a_1002 conda-forge xorg-renderproto 0.11.1 h14c3975_1002 conda-forge xorg-xextproto 7.3.0 h14c3975_1002 conda-forge xorg-xproto 7.0.31 h14c3975_1007 conda-forge xz 5.2.5 h516909a_1 conda-forge yaml 0.2.5 h516909a_0 conda-forge zarr 2.4.0 py_0 conda-forge zict 2.0.0 py_0 conda-forge zipp 3.1.0 py_0 conda-forge zlib 1.2.11 h516909a_1007 conda-forge zstd 1.4.5 h6597ccf_2 conda-forge ``` ### pip list ``` ```
greglucas commented 4 years ago

Yes, Matplotlib is deprecating flat shading because it was silently dropping the final row/column of data when the shapes were equal. With the update, shading="auto" does some magic under the hood when shape(X) = shape(Y) = shape(C) to calculate differences and increase X/Y by midpoints between the coordinates.

See here for a good description: https://matplotlib.org/3.3.1/gallery/images_contours_and_fields/pcolormesh_grids.html

Unfortunately, this dataset is not strictly increasing and so the difference between coordinates jumps from 360 to 0 in the middle of the xc array: np.diff(rasm.xc).max() # 359.8 which produces non-midpoint x-coordinates and auto is not doing what we'd expect.

You can take Cartopy away completely and just plot in x/y coordinates and the problem is evident there too. Figure_1

mathause commented 4 years ago

Then this is an upstream issue - thanks. Now the quest begins to find out if there is already a mpl issue on this. Feel free to close.

I add an example that works without xarray:

import numpy as np
import matplotlib.pyplot as plt
lon = np.concatenate((np.arange(180, 360, 5), np.arange(0, 180, 5)))
lat = np.arange(90, 0, -1)
data = np.random.randn(*lat.shape + lon.shape) + lon

# set some data to nan to make the problem more visible
data[::5, 35] = np.NaN

plt.pcolormesh(lon, lat, data, shading="auto")
mathause commented 4 years ago

Ahh but in pure matplotlib this does not work with shading="flat" either.

rasm.isel(time=1).Tair.plot.pcolormesh.plot(x="xc", y="yc", shading="flat")

rasm.isel(time=1).Tair.plot.pcolormesh.plot(x="xc", y="yc", shading="nearest")

So mpl can probably assume that the coordinates are monotonic.

greglucas commented 4 years ago

pcolormesh is defined to draw quadrilaterals where you specify the edges, not the midpoints. In addition, it assumes the end of one cell is the beginning of the next. See the diagram in the documentation https://matplotlib.org/3.3.1/api/_as_gen/matplotlib.axes.Axes.pcolormesh.html#matplotlib.axes.Axes.pcolormesh which I think does a pretty good job of describing that.

So, the coordinates don't have to be monotonic, but the cells will all be "connected" in a sense. If you open up an issue, feel free to ping me on it.

I think that a warning should probably be added to Matplotlib for the new 'auto' case when things are not strictly increasing/decreasing in the coordinate arrays and they try to find midpoints.

jklymak commented 4 years ago

I closed the upstream issue just so we can hash this out here. Of course Matplotlib is willing to make changes, particularly if the new default is going to end up being problematic.

I'm having trouble seeing what the way forward here is. The dataset above has a wrapping issue with no Cartopy projection in either "flat" or "nearest". Looking at the data set, that seems to be because the wrap occurs for different rows in different columns of xc and yc. That makes the cartesian quadrilaterals at that column wrap and leads to the distortion?

I assume that flat works with NorthPolarStereo because the data set was developed for that projection.

If we add a warning in Matplotlib, is it really just for "nearest" that there is an issue? I guess the "bad" thing we do for nearest is add a data point, which in this case is incorrect. I'm wondering if we should really do the grid re-centring in axes space and then transform back to data space...

greglucas commented 4 years ago

The thing I am still confused by is how the original xarray example ever worked in the previous versions because the wrap point was in the middle of the array and not one the values that was clipped. I'm curious if Cartopy previously did just fine transforming the coordinates and then identified those cells as wrapped and masked them out. So, Cartopy was essentially hiding the previous pcolormesh behavior.

The original example in pure MPL still is not correct even with MPL 3.2

import matplotlib.pyplot as plt
import xarray as xr
rasm = xr.tutorial.load_dataset("rasm")

ax = plt.subplot(111)

rasm.isel(time=1).Tair.plot.pcolormesh(
    ax=ax, x="xc", y="yc", shading='flat')

plt.show()
jklymak commented 4 years ago

The grid in that example looks like this in Cartesian co-ordinates (every 10th row and column)

image

If you form a quadrilateral in cartesian co-ordinates (default no cartopy) you have the ones near the edge with at least one corner that wraps. Presumably when you plot these from `NorthPolarStereo' there is no wrap.

fig, ax = plt.subplots()
ax.plot(rasm.xc[::10, ::10], rasm.yc[::10, ::10], '.', markersize=2);
ax.plot(rasm.xc[::10, 10], rasm.yc[::10, 10], markersize=2);
ax.plot(rasm.xc[::10, 230], rasm.yc[::10, 230], markersize=2);
ax.plot(rasm.xc[130, 230], rasm.yc[130, 230], 'd')
ax.plot(rasm.xc[140, 230], rasm.yc[140, 230], 'd')
ax.plot(rasm.xc[140, 240], rasm.yc[140, 240], 'd')
ax.plot(rasm.xc[130, 240], rasm.yc[130, 240], 'd')
jklymak commented 4 years ago

Just as prior art, xarray has _infer_interval_breaks which bails if the data is non-monotonic (unless over-ridden). I guess matplotlib could just do that. but - it would be nice to do the right thing if possible.

greglucas commented 4 years ago

I gave this a bit of thought and am not sure where we should fix this, either the MPL or Cartopy side.

Cartopy is currently using a private method from MPL _pcolorargs, that is pretty powerful and does a lot of error/bounds checking for us. I'm not sure we want to reproduce all of that in Cartopy directly because then we could just get out of sync in the future pretty easily with a new update. Maybe there is a way we could push some of that logic into a cbook.make_grids() public function or something similar?

I think the easiest way out is to add a keyword to _pcolorargs(..., wrap=None). The keyword would only be used by Cartopy. We could then account for the angle wrap inside of the grid generation routine directly. Change this: https://github.com/matplotlib/matplotlib/blob/682836cd084499ec9631098635adef26ba5d999e/lib/matplotlib/axes/_axes.py#L5684 to:

                        dX = np.diff(X, axis=1)
                        if wrap is not None:
                            # Account for angle wraps
                            dX = (dX + wrap/2) % wrap - wrap/2
                        dX = dX/2.

Any other ideas on how we can keep Cartopy and MPL in sync here?

jklymak commented 4 years ago

@greglucas There are two types of wraps here, if I understand correctly.

The simplest is what happens at -179 to +179 where there is mathematically a wrap, but after applying the projection there is no wrap i.e. in screen space the two data points can be next to each other (depending on central longitude of course). I think matplotlib could handle that wrap by doing the grid re-interpolation (or recentering) in screen space (i.e. after the projection transform has been applied) rather than before, which it does now. This would subtly skew existing grids for non-linear projections, but I don't think it would be too bad.

Your proposed wrap kwarg seems possible as a another solution to this - we'd need wrapX and wrapY. I guess I prefer doing the re-interpolation in screen space as thats more general. However, it could be more expensive!

The second type of wrapping is due to the user not giving you the data properly sorted into physical space. eg if the user gives you data from -180 to +180, but then makes the central longitude 80 degrees. The axes now go from -100 to -100, and data at -99 is not contiguous on screen with -101, but it is in the data set. In that case I think cartopy is on its own because matplotlib can't re-order the data or insert masks to create gaps.

greglucas commented 4 years ago

Correct, I think we are on the same page with multiple wraps.

Cartopy handles the "screen"/"projection" wraps currently. It can handle the points [181, 360] and map those to the negatives [-179, 0] and masks the the values that cross the wrap boundary and internally redraws those quads.

The part that is messing Cartopy up currently is handling inserting the new points to make the quad edges that aren't correct in the data-space. So, really all we need to do is add in proper data-space points and then Cartopy will handle the rest. That is the easy fix and will work with data going [179, 181, 183] and [359, 1, 3].

The harder case is actually if a user decides to put their data in this order: [179, -179, -177]. Cartopy even previously fails with that because it will draw the full quad just in the backwards orientation going from 179 to -179 covering most of the map and there is no boundary intersection/wrap.

I would have to give some more thought to how we could identify and mask quads if we do the interpolation in screen space instead. I am generally more of a proponent of doing the interpolations in the space that the data is given, not the one it is being mapped to. If a user gives us linearly spaced points, but wants to display them on a log-scaled axis do we add the bin edges in log-space (screen coordinates) or linear-space (data coordinates)? I'm not sure there is a correct answer here, really it probably boils down to preference.

jklymak commented 4 years ago

If a user gives us linearly spaced points, but wants to display them on a log-scaled axis do we add the bin edges in log-space (screen coordinates) or linear-space (data coordinates)? I'm not sure there is a correct answer here, really it probably boils down to preference.

I generally agree with you, but keep in mind its just a slightly different dx we are talking about - the quad centres will still be spaced as before, just their extents would be a bit different. To me that is better than completely breaking if a non-linear projection is used. If a user is fussy that the quad edges line up exactly 0.5 of the cell width in data space, then they can specify the quad edges explicitly and not rely on whatever we do.

greglucas commented 4 years ago

I just pushed up an option to handle this on Cartopy's side in #1646. That handles defining bin edges in data space before getting to Matplotlib. I think that would be agnostic of whether you change MPL to do screenspace interpolation or not. I guess that I find this way cleaner personally, but it does add some of the MPL grid interpolation code into Cartopy now.

jklymak commented 4 years ago

Without checking exactly what you did I think that's a great solution for downstream: don't count on matplotlib generic edge interpolation if you don't need to.

jklymak commented 4 years ago

OK, so for clarity Matplotlib has decided to not deal with fancy re-interpolation. If users encounter non-monotonic grid, indicative of a possible wrap issue, they now get a warning and are asked to specify the bin edges rather than the centers. (Thanks https://github.com/matplotlib/matplotlib/pull/18398 @greglucas!). If we continue to get feedback that this is inadequate we may consider trying to do something fancier like interpolation in display space, but it wasn't clear if people want quadrilaterals to be square in display space of data space, so we decided to punt for matplotlib 3.3.2.

greglucas commented 4 years ago

Sounds good. I think the new warning message is pretty clear now. And hopefully it is only a small number of people that will ever hit that path in the first place...

jklymak commented 4 years ago

Hmm, I think a lot of people will hit it. I plot pcolormesh all the time that double back on themselves. e.g. we have an underwater glider trying to drive west, but its getting pushed around some by currents, so sometimes goes backwards. A pcolormesh is perfectly reasonable in this case, and its perfectly reasonable to not have the x co-ordinate be monotonic.

greglucas commented 4 years ago

Interesting use case, neat. I generally bin/grid my data so it always ends up in some sorted order. I guess you're just OK with overlapping quads and hiding the 'lower' one, assuming your most recent data is best?

jklymak commented 4 years ago

One hopes the data aren't that different, and of course you can always plot in time if you want to unfold it. Binning the data is OK for some applications, but involves interpolation and/or averaging which is not always a luxury the data can afford.