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

PlateCarree projection errors cause gridline and view limit problems when cental longitude != 0 #1751

Open
brenmous opened this issue Mar 19, 2021 · 3 comments

Comments

@brenmous
Copy link

Description

I've been working on a global map centered on Australia (longitude 132) by using a PlateCarree projection centered on 132 for the axes. I set the extent using a 0-centered PlateCarree CRS, using coordinates relative to it to center on 132 (west = -48, east = 312, south = -90, north = 90).

Gridline labels worked fine with a 0 centered CRS, but not with the 132 CRS. The right labels are missing (bottom labels missing is intentional):

cartopy_bug_no_labels

I found that point transforms when using PlateCarree not centered on 0 give strange results:

>>> cl = cp.crs.PlateCarree(central_longitude=0)
>>> cl.transform_point(180, 0, cl)
(180.0, 0.0)
>>> cl.transform_point(-180, 0, cl)
(-180.0, 0.0)
>>> cl132 = cp.crs.PlateCarree(central_longitude=132)
>>> cl132.transform_point(180, 0, cl132)
(-179.99999999999997, 0.0)
>>> cl132.transform_point(-180, 0, cl132)
(-180.0, 0.0)

180 will always transform to -180, along with some floating point error.

This line

line = self.axes.projection.transform_points(

in gridliner projects the latitude and longitude gridlines to the axes CRS. As I'm using the 132 centered CRS for both the axes and gridlines, it transforms against itself. This turns the start and end of the longitude lines to -180, -180, rather than -180, 180. So any labels intersecting 180 actually get placed at -180 (or at least I think they do, I expected setting left_labels to False and right_labels to True would cause the left labels to appear, but it doesn't - either way it stops the right labels from being drawn).

My workaround for this was to only perform the projection if the gridliner CRS != axes CRS. This gives me the right labels, but I'm still missing the rightmost longitude label:

cartopy_bug_no_last_lon_label

This is due to precision error when set_extent projects the extent values:

map_projection = cp.crs.PlateCarree(central_longitude=0)
print(f"View extent: l={l}, r={r}, b={b}, t={t}")
ax.set_extent([l, r, b, t], crs=map_projection)
print(f"viewLim: {ax.viewLim}")

# Outputs:
# View extent: l=-48.0, r=312.0, b=-90.0, t=90.0
# viewLim: Bbox(x0=-180.0, y0=-90.0, x1=179.999656677246, y1=90.0)

In this case there's precision error transforming between the 0 CRS and the 132 CRS.

The longitude gridlines now correctly span -180 to 180, but because the right-hand map boundary is less than 180, the longitude gridline doesn't intersect with the map

if line.intersects(map_boundary):

so the label doesn't get drawn. The solution to this has been to set xlim and ylim manually (set_global also works provided you want 90 south to 90 north):

cartopy_bug_desired_result

With that said, thanks for working on Cartopy. Besides the odd issue here and there, it's been really great to use.

Code to reproduce

import cartopy as cp
import matplotlib.pyplot as plt

cl132 = cp.crs.PlateCarree(central_longitude=132)
cl0 = cp.crs.PlateCarree(central_longitude=0)

# Global extents for a map centred on lon 132 relative to central lon 0
west, south, east, north = -48, -90, 312, 90

ax = plt.axes(projection=cl132)
ax.set_extent((west, east, south, north), crs=cl0)
ax.add_feature(cp.feature.LAND)
ax.gridlines(draw_labels=True, crs=cl132)

plt.show()

Operating system

Ubuntu 18.04

Cartopy version

0.18, Python 3.7.9

pip list

affine                        2.3.0
alabaster                     0.7.12
altgraph                      0.17
attrs                         20.3.0
awscli                        1.18.166
Babel                         2.8.0
botocore                      1.19.6
Cartopy                       0.18.0
certifi                       2020.11.8
chardet                       3.0.4
click                         7.1.2
click-plugins                 1.1.1
cligj                         0.7.0
colorama                      0.4.3
cx-Freeze                     6.5.3
cycler                        0.10.0
declxml                       1.1.3
decorator                     4.4.2
dill                          0.3.3
docutils                      0.16
eatws-skip-client             0.12.9
Flask                         1.1.2
Flask-Cors                    3.0.9
Flask-WTF                     0.14.3
future                        0.18.2
gprof2dot                     2019.11.30
gunicorn                      20.0.4
idna                          2.10
imagesize                     1.2.0
importlib-metadata            2.0.0
iniconfig                     1.1.1
itsdangerous                  1.1.0
Jinja2                        2.11.2
jmespath                      0.10.0
kiwisolver                    1.3.1
lxml                          4.6.1
MarkupSafe                    1.1.1
matplotlib                    3.3.2
mock                          4.0.2
mypy                          0.790
mypy-extensions               0.4.3
networkx                      2.5
Nuitka                        0.6.12.1
numpy                         1.19.4
obspy                         1.2.2
packaging                     20.4
Pillow                        8.0.1
pip                           20.2.3
pip-tools                     5.3.1
pluggy                        0.13.1
py                            1.9.0
pyasn1                        0.4.8
Pygments                      2.7.2
pyinstaller                   4.2
pyinstaller-hooks-contrib     2020.11
pyparsing                     2.4.7
pyprof2calltree               1.4.5
PyQt5                         5.15.1
PyQt5-sip                     12.8.1
pyshp                         2.1.2
pytest                        6.1.2
python-dateutil               2.8.1
pytz                          2020.4
PyYAML                        5.3.1
rasterio                      1.1.8
requests                      2.24.0
requests-toolbelt             0.9.1
rsa                           4.5
s3transfer                    0.3.3
schedule                      0.6.0
scipy                         1.5.4
setuptools                    54.1.2
Shapely                       1.7.1
six                           1.15.0
snowballstemmer               2.0.0
snuggs                        1.4.7
Sphinx                        3.3.0
sphinx-autodoc-typehints      1.11.1
sphinx-click                  2.6.0
sphinx-rtd-theme              0.5.1
sphinxcontrib-applehelp       1.0.2
sphinxcontrib-devhelp         1.0.2
sphinxcontrib-htmlhelp        1.0.3
sphinxcontrib-jsmath          1.0.1
sphinxcontrib-programoutput   0.16
SQLAlchemy                    1.3.20
toml                          0.10.2
tqdm                          4.57.0
typed-ast                     1.4.1
typing                        3.7.4.3
typing-extensions             3.7.4.3
urllib3                       1.25.11
Werkzeug                      1.0.1
WTForms                       2.3.3
xmltodict                     0.12.0
zipp                          3.4.0
brenmous added a commit to brenmous/cartopy that referenced this issue Mar 21, 2021
brenmous added a commit to brenmous/cartopy that referenced this issue Mar 21, 2021
@brenmous
Copy link
Author

Have tested on the master branch and still exists, also upgrading Proj to 8.0.0 didn't fix it (was using 5.1.0 before).

@dopplershift
Copy link
Contributor

@brenmous I'm not sure how you're building and compiling CartoPy, but I'm 99% sure that compiling CartoPy against Proj 8.0 is impossible. (See #1289 #1742 #1752)

@brenmous
Copy link
Author

brenmous commented Mar 22, 2021

You're correct, I hadn't reinstalled Cartopy properly - some pj_get_version function or similar can't be found with Proj 8.0.

@dopplershift so having rebuilt Cartopy against 7.2.1 (and checking _crs.PROJ4_RELEASE to make sure I did it properly):

  • branch v0.18.0 + Proj 5.2: error present
  • branch v0.18.0 + Proj 7.2.1: error present
  • master + Proj 5.2: right-hand latitude labels appear but right-most longitude labels still missing
  • master + Proj 7.2.1: same as Proj 5.2

Code on the master branch fixes the missing latitude labels. The error in projecting the view limit still exists (x1=179.999656677246 across all configurations I tried).

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

No branches or pull requests

2 participants