cartopy
cartopy copied to clipboard
SciTools
The bounds given for a Geostationary projection do not plot correctly
Description
This is a bug with the plot range that is shown in cartopy. With a Geostationary projection if you use the bounds given in the projection_object.boundary.bounds to set your extent then the plot will show only a small rectangle of the middle of the globe. If you remove a very small amount from each side of the bounds, suddenly it shows the whole globe again.
There is no crash or traceback message, the program simply creates incorrect plots.
Code to reproduce
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import cartopy.crs as crs
import cartopy
# This geostationary projection matches that of GOES 16.
proj_object = crs.Geostationary(central_longitude=-75.0, satellite_height=35786023.0, sweep_axis='x', )
bounding_axes = proj_object.boundary.bounds
# note, the bounds are in the wrong order, so reorder them
bounding_axes = (bounding_axes[0], bounding_axes[2], bounding_axes[1], bounding_axes[3])
# plot a figure that shows the broken bounds behavior
fig = plt.figure( )
axes = fig.add_subplot(111, projection=proj_object)
axes.set_extent(bounding_axes, crs=proj_object)
axes.add_feature(cartopy.feature.LAND, facecolor="green", ) # not nessicary, but makes it clearer where we are
fig.savefig("./fdtest1.png", )
plt.close(fig)
# plot a figure that works around the bug
bounds_offset = 1.0
new_bounds = (bounding_axes[0] + bounds_offset, bounding_axes[1] - bounds_offset,
bounding_axes[2] + bounds_offset, bounding_axes[3] - bounds_offset,)
fig = plt.figure( )
axes = fig.add_subplot(111, projection=proj_object)
axes.set_extent(new_bounds, crs=proj_object)
axes.add_feature(cartopy.feature.LAND, facecolor="green", ) # not nessicary, but makes it clearer where we are
fig.savefig("./fdtest2.png", )
plt.close(fig)
Here are the plots I get:
Traceback
(No traceback is created.)
Full environment definition
Operating system
MacBook Pro with MacOS 12.6.7
Cartopy version
cartopy 0.21.1
conda list
% mamba list
# packages in environment at /Users/evas/mambaforge/envs/aitf_ql:
#
# Name Version Build Channel
aitf-ql 0.8 dev_0 <develop>
appdirs 1.4.4 pyhd3eb1b0_0
blas 1.0 openblas
brotli 1.0.9 hca72f7f_7
brotli-bin 1.0.9 hca72f7f_7
brotlipy 0.7.0 py38h9ed2024_1003
bzip2 1.0.8 h1de35cc_0
c-ares 1.19.0 h6c40b1e_0
ca-certificates 2023.05.30 hecd8cb5_0
cartopy 0.21.1 py38hfbf152b_0
certifi 2023.7.22 py38hecd8cb5_0
cffi 1.15.1 py38h6c40b1e_3
cftime 1.6.2 py38hacda100_0
charset-normalizer 2.0.4 pyhd3eb1b0_0
contourpy 1.0.5 py38haf03e11_0
cryptography 41.0.2 py38h3b477ad_0
cycler 0.11.0 pyhd3eb1b0_0
fonttools 4.25.0 pyhd3eb1b0_0
freetype 2.12.1 hd8bbffd_0
geos 3.8.0 hb1e8313_0
giflib 5.2.1 h6c40b1e_3
hdf4 4.2.13 h39711bb_2
hdf5 1.12.1 ha01d115_3
idna 3.4 py38hecd8cb5_0
importlib_resources 5.2.0 pyhd3eb1b0_1
jpeg 9e h6c40b1e_1
kiwisolver 1.4.4 py38hcec6c5f_0
krb5 1.20.1 h428f121_1
lcms2 2.12 hf1fd2bf_0
lerc 3.0 he9d5cce_0
libbrotlicommon 1.0.9 hca72f7f_7
libbrotlidec 1.0.9 hca72f7f_7
libbrotlienc 1.0.9 hca72f7f_7
libcurl 8.1.1 hf20ceda_2
libcxx 14.0.6 h9765a3e_0
libdeflate 1.17 hb664fd8_0
libedit 3.1.20221030 h6c40b1e_0
libev 4.33 h9ed2024_1
libffi 3.4.4 hecd8cb5_0
libgfortran 5.0.0 11_3_0_hecd8cb5_28
libgfortran5 11.3.0 h9dfd629_28
libnetcdf 4.8.1 h9f5a9a2_4
libnghttp2 1.52.0 h9beae6a_1
libopenblas 0.3.21 h54e7dc3_0
libpng 1.6.39 h6c40b1e_0
libssh2 1.10.0 h04015c4_2
libtiff 4.5.0 hcec6c5f_2
libwebp 1.2.4 hf6ce154_1
libwebp-base 1.2.4 h6c40b1e_1
libzip 1.8.0 h29ab7a1_1
llvm-openmp 14.0.6 h0dcd299_0
lz4-c 1.9.4 hcec6c5f_0
matplotlib-base 3.7.1 py38hda11e5a_1
munkres 1.1.4 py_0
ncurses 6.4 hcec6c5f_0
netcdf4 1.6.2 py38hd243f81_0
numpy 1.24.3 py38h57a7bef_0
numpy-base 1.24.3 py38hc93c6d9_0
openssl 3.0.10 hca72f7f_0
packaging 23.0 py38hecd8cb5_0
pillow 9.4.0 py38hcec6c5f_0
pip 23.2.1 py38hecd8cb5_0
pooch 1.4.0 pyhd3eb1b0_0
proj 8.2.1 hd69def0_0
pycparser 2.21 pyhd3eb1b0_0
pyopenssl 23.2.0 py38hecd8cb5_0
pyparsing 3.0.9 py38hecd8cb5_0
pyproj 3.4.1 py38hf797ea4_0
pyshp 2.1.3 pyhd3eb1b0_0
pysocks 1.7.1 py38_1
python 3.8.17 h5ee71fb_0
python-dateutil 2.8.2 pyhd3eb1b0_0
readline 8.2 hca72f7f_0
requests 2.31.0 py38hecd8cb5_0
scipy 1.10.1 py38h9034365_1
setuptools 68.0.0 py38hecd8cb5_0
shapely 2.0.1 py38h797b0b3_0
six 1.16.0 pyhd3eb1b0_1
sqlite 3.41.2 h6c40b1e_0
tk 8.6.12 h5d9f67b_0
urllib3 1.26.16 py38hecd8cb5_0
wheel 0.38.4 py38hecd8cb5_0
xz 5.4.2 h6c40b1e_0
zipp 3.11.0 py38hecd8cb5_0
zlib 1.2.13 h4dc903c_0
zstd 1.5.5 hc035e20_0
pip list
% pip list
Package Version Editable project location
------------------- --------- ---------------------------------------------------------
aitf-ql 0.8 /Users/evas/Dev/CSPP/CSPP FW Quicklooks/AIT_FW_Quicklooks
appdirs 1.4.4
brotlipy 0.7.0
Cartopy 0.21.1
certifi 2023.7.22
cffi 1.15.1
cftime 1.6.2
charset-normalizer 2.0.4
contourpy 1.0.5
cryptography 41.0.2
cycler 0.11.0
fonttools 4.25.0
idna 3.4
importlib-resources 5.2.0
kiwisolver 1.4.4
matplotlib 3.7.1
munkres 1.1.4
netCDF4 1.6.2
numpy 1.24.3
packaging 23.0
Pillow 9.4.0
pip 23.2.1
pooch 1.4.0
pycparser 2.21
pyOpenSSL 23.2.0
pyparsing 3.0.9
pyproj 3.4.1
pyshp 2.1.3
PySocks 1.7.1
python-dateutil 2.8.2
requests 2.31.0
scipy 1.10.1
setuptools 68.0.0
shapely 2.0.1
six 1.16.0
urllib3 1.26.16
wheel 0.38.4
zipp 3.11.0
I think there's something wrong (maybe some precision issue) with the interpolator going on here. If i increase the number of points in the boundary for Geostationary (to say 1001 from 91), I can actually trigger a ValueError: Axis limits cannot be NaN or Inf. The points in that boundary itself, though, all seem to be valid since (using pyproj.Transform directly) they all transform fine to PlateCarre/Geodetic--no nans. Not sure why adding points produced worse results here.
In the meantime, there are a couple easy enough work-arounds. If all you want is to just have the limits show the entire valid range of the projection, just use axes.set_global(). Alternatively, since you directly have x/y ranges in projected space, you can use set_xlim/set_ylim with those:
x1, x2, y1, y2 = bounding_axes
axes.set_xlim(x1, x2)
axes.set_ylim(y1, y2)
These both avoid the problematic (and in this case, completely unnecessary) coordinate projection/transformation pipeline.
I would bet a substantial amount of money that this is because the extent is represented as a LineString and what the interpolator needs to sample over is a Polygon:
https://github.com/SciTools/cartopy/blob/d1b43fbf641729e68b5a31da38d34d06bd06b593/lib/cartopy/mpl/geoaxes.py#L833
This issue is something I addressed in #1739 .
