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 365 forks source link

Differences between when transform_first is True and False #2233

Open malloryprow opened 1 year ago

malloryprow commented 1 year ago

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 transform_first_False transform_first_True

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 ```
lgolston commented 1 year ago

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?

malloryprow commented 1 year ago

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.

rcomer commented 1 year ago

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

lgolston commented 1 year ago

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: transform_first_True_no_sort

malloryprow commented 1 year ago

Definitely much better. Maybe this is just a bad grid transformation to do for using transform_first=True?

lgolston commented 1 year ago

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. transform_first_True_no_sort_add_nan

malloryprow commented 1 year ago

Thanks @lgolston!!

greglucas commented 1 year ago

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.