Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

The bounds given for a Geostationary projection do not plot correctly #2240

Open
evas-ssec opened this issue Aug 29, 2023 · 2 comments
Open

Comments

@evas-ssec
Copy link

evas-ssec commented Aug 29, 2023

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

fdtest2

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
@dopplershift
Copy link
Contributor

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.

@kdpenner
Copy link
Contributor

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:

domain_in_crs = sgeom.polygon.LineString([[x1, y1], [x2, y1],

This issue is something I addressed in #1739 .

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

No branches or pull requests

3 participants