cartopy
cartopy copied to clipboard
SciTools
Differences between when transform_first is True and False
Description
I'm using Cartopy to create spatial maps of 24 hour precipitation accumulations from NCEP's GFS model. I'm plotting the data that is on a 0.25 degree lat-lon grid (cartopy.crs.PlateCarree) to Lambert Conformal (cartopy.crs.LambertConformal) projection. I started noticing some oddities in the images. After some testing, I traced it back to the value of transform_first. With transform_first=False things look fine; with transform_first=True it looks like something odd is happening. I would think the images would look the same as I thought transform_first was a setting to help speedup runtime. This appears to be a bug in transform_first=True.
Attached are two images displaying the issue. Additionally text files with the data needed to use the code.
gfs_init2023080412_f24_APCP_A24.txt
gfs_init2023080412_f24_lat.txt
gfs_init2023080412_f24_lon.txt
Code to reproduce
#!/usr/bin/env python3
'''
Name: global_det_atmos_plots_precip_spatial_map.py
Contact(s): Mallory Row
Abstract: This script generates a spatial map for precip
'''
import netCDF4 as netcdf
import pyproj
import numpy as np
import matplotlib
matplotlib.use('agg')
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.util import add_cyclic_point
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
from cartopy import config
# Set contour levels and color map
clevs = [0.01, 0.10, 0.25, 0.50, 0.75, 1.00, 1.25, 1.50,
1.75, 2.00, 2.50, 3.00, 4.00, 5.00, 7.00, 10.00,
15.00, 20.00]
colorlist = ['#7fff00', '#00cd00', '#008b00', '#104e8b',
'#1e90ff', '#00b2ee', '#00eeee', '#8968cd',
'#912cee', '#8b008b', '#8b0000', '#cd0000',
'#ee4000', '#ff7f00', '#cd8500', '#ffd700',
'#ffff00', '#ffff02']
cmap_over_color = '#ffaeb9'
cmap = matplotlib.colors.ListedColormap(colorlist)
cmap_over_color = cmap_over_color
norm = matplotlib.colors.BoundaryNorm(clevs, cmap.N)
# Read in data
precip_lat = np.loadtxt('gfs_init2023080412_f24_lat.txt')
precip_lon = np.loadtxt('gfs_init2023080412_f24_lon.txt')
precip_APCP_A24 = np.loadtxt('gfs_init2023080412_f24_APCP_A24.txt')
x, y = np.meshgrid(precip_lon, precip_lat)
precip_APCP_A24 = precip_APCP_A24 * 0.0393701 # convert to inches
# Set projection
myproj=ccrs.LambertConformal()
# Make transform_first=False
fig = plt.figure(figsize=(8.,6.))
gs_hspace, gs_wspace = 0, 0
gs_bottom, gs_top = 0.125, 0.85
gs = gridspec.GridSpec(1,1, bottom=gs_bottom, top=gs_top,
hspace=gs_hspace, wspace=gs_wspace)
ax1 = fig.add_subplot(gs[0], projection=myproj)
ax1.add_feature(cfeature.COASTLINE.with_scale('50m'),
zorder=2, linewidth=1)
ax1.add_feature(cfeature.BORDERS.with_scale('50m'),
zorder=2, linewidth=1)
ax1.add_feature(cfeature.STATES.with_scale('50m'),
zorder=2, linewidth=1)
fig.suptitle('transform_first=False')
CF1 = ax1.contourf(x, y, precip_APCP_A24,
transform=ccrs.PlateCarree(),
levels=clevs, norm=norm,
cmap=cmap, extend='max',
transform_first=False)
CF1.cmap.set_over(cmap_over_color)
plt.savefig('transform_first_False.png')
# Make transform_first=True
fig = plt.figure(figsize=(8.,6.))
gs_hspace, gs_wspace = 0, 0
gs_bottom, gs_top = 0.125, 0.85
gs = gridspec.GridSpec(1,1, bottom=gs_bottom, top=gs_top,
hspace=gs_hspace, wspace=gs_wspace)
ax1 = fig.add_subplot(gs[0], projection=myproj)
ax1.add_feature(cfeature.COASTLINE.with_scale('50m'),
zorder=2, linewidth=1)
ax1.add_feature(cfeature.BORDERS.with_scale('50m'),
zorder=2, linewidth=1)
ax1.add_feature(cfeature.STATES.with_scale('50m'),
zorder=2, linewidth=1)
fig.suptitle('transform_first=True')
CF1 = ax1.contourf(x, y, precip_APCP_A24,
transform=ccrs.PlateCarree(),
levels=clevs, norm=norm,
cmap=cmap, extend='max',
transform_first=True)
CF1.cmap.set_over(cmap_over_color)
plt.savefig('transform_first_True.png')
plt.close('all')
Traceback
Full environment definition
Operating system
NCEP's WCOSS2
Cartopy version
0.21.1
conda list
pip list
Package Version
------------------ ----------
asttokens 2.2.1
attrs 22.1.0
backcall 0.2.0
Cartopy 0.21.1
certifi 2022.9.24
cftime 1.6.2
charset-normalizer 3.1.0
click 8.1.3
contourpy 1.0.5
cycler 0.11.0
dateutils 0.6.12
decorator 5.1.1
distlib 0.3.6
executing 1.2.0
filelock 3.8.0
Flask 2.2.3
fonttools 4.37.4
geos 0.2.3
idna 3.4
imageio 2.22.1
iniconfig 1.1.1
ipython 8.9.0
itsdangerous 2.1.2
jedi 0.18.2
Jinja2 3.1.2
joblib 1.2.0
kiwisolver 1.4.4
lxml 4.9.2
MarkupSafe 2.1.2
matplotlib 3.7.1
matplotlib-inline 0.1.6
metplus 4.1.4
MetPy 1.4.1
netCDF4 1.6.3
numpy 1.25.2
packaging 21.3
pandas 1.5.3
parso 0.8.3
pexpect 4.8.0
pickleshare 0.7.5
Pillow 9.2.0
Pint 0.20.1
pip 23.1.2
pipdeptree 2.3.1
pipenv 2022.10.11
platformdirs 2.5.2
plotly 5.14.1
pluggy 1.0.0
pooch 1.7.0
prompt-toolkit 3.0.36
ptyprocess 0.7.0
pure-eval 0.2.2
py 1.11.0
Pygments 2.14.0
pyparsing 3.0.9
pyproj 3.4.0
pyshp 2.3.1
pytest 7.1.3
python-dateutil 2.8.2
pytz 2022.4
PyYAML 6.0
requests 2.28.2
scikit-learn 1.2.2
scipy 1.10.1
setuptools 68.0.0
shapely 2.0.1
six 1.16.0
stack-data 0.6.2
tenacity 8.2.2
threadpoolctl 3.1.0
tomli 2.0.1
traitlets 5.9.0
tzdata 2023.3
urllib3 1.26.15
virtualenv 20.16.5
virtualenv-clone 0.5.7
wcwidth 0.2.6
Werkzeug 2.2.3
wheel 0.40.0
xarray 2023.3.0
Unfortunately I think some loss of accuracy associated with the speedup from transform_first=True is expected - this is noted at the contour docs and associated example, although the differences are pretty significant in this case. Maybe the warning in the docs should be strengthened / better clarify the appropriate use case?
That might be nice. A warning perhaps if using transform_first=True? Yes like you said the differences are quite significant in this case for whatever reason.
The documentation for contourf says that both x and y must be monotonic. There is logic in the tranform_first handling to ensure that x is monotonic, but I wonder if it's actually possible for both to be with the current projection?
https://github.com/SciTools/cartopy/blob/2e722fc93221a60787ec31af80f0b32e637ae0d0/lib/cartopy/mpl/geoaxes.py#L352-L355
Great point @rcomer! I have investigated a little bit and that is definitely causing a problem - at low and high projected (and sorted) x values, y is jumping back and forth causing the vertical striping. There are still artifacts (notably, the continents still get shaded green), but simply disabling the sort otherwise creates a more sensible output here:
Definitely much better. Maybe this is just a bad grid transformation to do for using transform_first=True?
I don't think there is an inherent reason why the fast version couldn't perform better here. Looking again, with sorting turned off there is a large discontinuity in x around index 336. If a column of nan are inserted there, the land/ocean problem is fixed. Now there only seems to be pixel level differences compared to the original. Will have to do more testing to see if this approach can generalize. I would still not recommend transform_first=true for a final version of a figure, but it is very convenient that it is running ~100x faster.
Perhaps lexsort applied to x and then y would be useful here instead of just argsorting the x values?
https://numpy.org/doc/stable/reference/generated/numpy.lexsort.html
However, the key difference between transform first and the default is that we are drawing the contours in projected coordinates. So with the projection having that wedge out of it at the top, the contours will "jump" across that gap and try to use values on either side of it which is likely not physical. There currently isn't a wrapping/jump handling contour algorithm in pycontour (that would be a great addition if someone is interested in adding it!) so there isn't a great way of handling these projections that have wedges out of them due to the underlying algorithms.
