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.43k stars 363 forks source link

pcolormesh in 0.18.0 takes about 10x time than in 0.17.0 #1676

Closed fercook closed 1 year ago

fercook commented 3 years ago

Description

I am creating a simple heatmap with pcolormesh from a netCDF file (*)

Using a new environment with Anaconda ( conda install matplotlib; conda install cartopy) installs matplotlib==3.3.2 and Cartopy==0.18.0. The code introduced below takes between 20 and 50 seconds to execute, depending on the machine (tested on 4, 2 Ubuntus and 2 Macs, different processors/ram/etc).

Using a new environment with conda install matplotlib==3.2.2; conda install cartopy==0.17.0, the code completes in between 1.5 and 3 seconds.

Just to make sure it's clear: Both versions work ok, but the newest one takes about 10x times to do it. Notice that to install 0.17.0 I had to downgrade matplotlib as there is an error with a function parameter. This means that I don't know if the performance issue is Cartopy or Matplotlib related.

(*) dataset file is available at https://github.com/fercook/ncexamplefile/blob/main/ERA5_REANALYSIS_air_temperature.nc, extracted from from Copernicus Data Store with the example script https://cds.climate.copernicus.eu/toolbox-editor/5196/21-calculate-regional-mean-and-anomalies

Code to reproduce

import matplotlib.pyplot as plt

from cartopy import config
import cartopy.crs as ccrs

from netCDF4 import Dataset as netcdf_dataset

from datetime import datetime

fname = "./ERA5_REANALYSIS_air_temperature.nc"
dataset = netcdf_dataset(fname)
tas = dataset.variables['tas'][:, :]
lats = dataset.variables['lat'][:]
lons = dataset.variables['lon'][:]

s=datetime.now()
ax = plt.axes(projection=ccrs.Robinson())
ax.pcolormesh(lons, lats, tas, transform = ccrs.PlateCarree())
plt.savefig('example.png', bbox_inches='tight')
e=datetime.now()
print("Time:",e-s)
Full environment definition ### Operating system Linux and macOS ### Cartopy version 0.18.0 vs 0.17.0 ### conda list ``` # packages in environment at /media/HDD/anaconda3hospitales/envs/mapping: # # Name Version Build Channel _libgcc_mutex 0.1 main blas 1.0 mkl brotlipy 0.7.0 py38h7b6447c_1000 bzip2 1.0.8 h7b6447c_0 ca-certificates 2020.10.14 0 cartopy 0.17.0 py38h3e2e9bd_1 certifi 2020.6.20 pyhd3eb1b0_3 cffi 1.14.3 py38he30daa8_0 cftime 1.2.1 py38heb32a55_0 chardet 3.0.4 py38_1003 cryptography 3.1.1 py38h1ba5d50_0 curl 7.71.1 hbc83047_1 cycler 0.10.0 py38_0 dbus 1.13.18 hb2f20db_0 expat 2.2.10 he6710b0_2 fontconfig 2.13.0 h9420a91_0 freetype 2.10.4 h5ab3b9f_0 geos 3.8.0 he6710b0_0 glib 2.66.1 h92f7085_0 gst-plugins-base 1.14.0 hbbd80ab_1 gstreamer 1.14.0 hb31296c_0 hdf4 4.2.13 h3ca952b_2 hdf5 1.10.4 hb1b8bf9_0 icu 58.2 he6710b0_3 idna 2.10 py_0 intel-openmp 2020.2 254 jpeg 9b h024ee3a_2 kiwisolver 1.3.0 py38h2531618_0 krb5 1.18.2 h173b8e3_0 lcms2 2.11 h396b838_0 ld_impl_linux-64 2.33.1 h53a641e_7 libcurl 7.71.1 h20c2e04_1 libedit 3.1.20191231 h14c3975_1 libffi 3.3 he6710b0_2 libgcc-ng 9.1.0 hdf63c60_0 libgfortran-ng 7.3.0 hdf63c60_0 libnetcdf 4.7.3 hb80b6cc_0 libpng 1.6.37 hbc83047_0 libssh2 1.9.0 h1ba5d50_1 libstdcxx-ng 9.1.0 hdf63c60_0 libtiff 4.1.0 h2733197_1 libuuid 1.0.3 h1bed415_2 libxcb 1.14 h7b6447c_0 libxml2 2.9.10 hb55368b_3 lz4-c 1.9.2 heb0550a_3 matplotlib 3.2.2 0 matplotlib-base 3.2.2 py38hef1b27d_0 mkl 2020.2 256 mkl-service 2.3.0 py38he904b0f_0 mkl_fft 1.2.0 py38h23d657b_0 mkl_random 1.1.1 py38h0573a6f_0 ncurses 6.2 he6710b0_1 netcdf4 1.5.3 py38hbf33ddf_0 numpy 1.19.2 py38h54aff64_0 numpy-base 1.19.2 py38hfa32c7d_0 olefile 0.46 py_0 openssl 1.1.1h h7b6447c_0 owslib 0.20.0 py_0 pcre 8.44 he6710b0_0 pillow 8.0.1 py38he98fc37_0 pip 20.2.4 py38_0 proj 7.0.1 h59a7b90_1 pycparser 2.20 py_2 pyepsg 0.4.0 py_0 pykdtree 1.3.1 py38hdd07704_1002 pyopenssl 19.1.0 py_1 pyparsing 2.4.7 py_0 pyproj 2.6.1.post1 py38h61f852b_1 pyqt 5.9.2 py38h05f1152_4 pyshp 2.1.2 py_0 pysocks 1.7.1 py38_0 python 3.8.5 h7579374_1 python-dateutil 2.8.1 py_0 pytz 2020.1 py_0 pyyaml 5.3.1 py38h7b6447c_1 qt 5.9.7 h5867ecd_1 readline 8.0 h7b6447c_0 requests 2.24.0 py_0 scipy 1.5.2 py38h0b6359f_0 setuptools 50.3.0 py38hb0f4dca_1 shapely 1.7.1 py38h98ec03d_0 sip 4.19.13 py38he6710b0_0 six 1.15.0 py_0 sqlite 3.33.0 h62c20be_0 tk 8.6.10 hbc83047_0 tornado 6.0.4 py38h7b6447c_1 urllib3 1.25.11 py_0 wheel 0.35.1 py_0 xz 5.2.5 h7b6447c_0 yaml 0.2.5 h7b6447c_0 zlib 1.2.11 h7b6447c_3 zstd 1.4.5 h9ceee32_0 ``` ### pip list ``` Package Version ------------------- ----------- adium-theme-ubuntu 0.3.4 appdirs 1.4.4 ccsm 0.9.12.3 certifi 2019.3.9 chardet 3.0.4 compizconfig-python 0.9.12.3 configparser 4.0.2 contextlib2 0.6.0.post1 cryptography 1.2.3 distlib 0.3.0 enum34 1.1.2 filelock 3.0.12 gyp 0.1 idna 2.8 importlib-metadata 1.6.0 importlib-resources 1.5.0 iotop 0.6 ipaddress 1.0.16 pathlib2 2.3.5 pip 20.0.2 pipenv 2018.11.26 pyasn1 0.1.9 pycrypto 2.6.1 pygobject 3.20.0 pyOpenSSL 0.15.1 requests 2.21.0 scandir 1.10.0 scour 0.32 setuptools 44.1.0 singledispatch 3.4.0.3 six 1.10.0 trash-cli 0.12.9.14 typing 3.7.4.1 unity-lens-photos 1.0 urllib3 1.24.2 virtualenv 20.0.20 virtualenv-clone 0.5.4 wheel 0.29.0 zipp 1.2.0 ```
greglucas commented 3 years ago

The root cause of this issue is that the shading in pcolormesh has changed from "flat" in MPL 3.2 to "auto" in 3.3, which will center your coordinates for you. https://matplotlib.org/3.3.0/gallery/images_contours_and_fields/pcolormesh_levels.html#centered-coordinates

The problem with Cartopy is that under the hood if you have point straddling the boundary, we will call pcolor to create those boundary cells. The pcolor calls are really slow. The simple way to speed this up is to do what the old Matplotlib did and truncate the data values tas[:-1, :-1] when calling pcolormesh.

Note that if you don't have any longitudes being wrapped, then the calls are just as fast.

import time

import matplotlib.pyplot as plt
import numpy as np
import cartopy.crs as ccrs

lats = np.linspace(-60, 60, 1000)
# -180, 180 is slow
# preventing wrapped coordinates is fast (-170, 170)
lons = np.linspace(-170, 170, 1000)
z = lats[:, np.newaxis] * lons[np.newaxis, :]

t0 = time.time()
fig = plt.figure()
ax = plt.axes(projection=ccrs.Robinson())
ax.pcolormesh(lons, lats, z, transform=ccrs.PlateCarree())
fig.canvas.draw()
t1 = time.time()
print("Time:", t1 - t0)
greglucas commented 1 year ago

Closing because I don't think there is anything that we can do within Cartopy for this. It is fundamentally Matplotlib expanding the grid and causing more path interpolation to take place.