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

Boundary splitting with a small Globe #1572

Open
kadrlica opened this issue May 21, 2020 · 1 comment
Open

Boundary splitting with a small Globe #1572

kadrlica opened this issue May 21, 2020 · 1 comment

Comments

@kadrlica
Copy link

kadrlica commented May 21, 2020

Description

Thanks for all the hard work on this nice tool.

I've started to transfer over some code that used basemap for plotting astronomical sky maps. Since I want a spherical Globe with no direct connection to physical scale of the Earth, I thought it would make sense to use a unit sphere:

globe = ccrs.Globe(ellipse=None,semimajor_axis=1.0,semiminor_axis=1.0)

This seems to work for the most part; however, when I try to draw a LineString that cross the global map boundary, it fails to break properly at the map boundary for some latitudes. I'm guessing that this is some numerical resolution issue, since this issues disappears as the radius of the globe is increased (and appears to be gone completely at radius = 1e5). My guess is that there is some implicit assumption of the globe size hidden somewhere, but I'm not sure where this can be found (or if I can somehow increase the resolution from the api directly).

Code to reproduce

import matplotlib.pyplot as plt
import numpy as np

import cartopy.crs as ccrs
from shapely.geometry.polygon import LineString

fig = plt.figure(figsize=(12,4))
plt.subplots_adjust(left=0.01,right=0.99,wspace=0.01)

radii = [1.0,1e4,1e5]
lon = np.linspace(0,360,100)
lat = 60*np.sin(np.radians(lon))

for i,radius in enumerate(radii):
    globe=ccrs.Globe(ellipse=None,semimajor_axis=radius,semiminor_axis=radius)
    projection=ccrs.Mollweide(globe=globe)
    geodetic = ccrs.Geodetic(globe=globe)

    ax = fig.add_subplot(1, len(radii), i+1, projection=projection)
    plt.title('radius = %g'%radius)

    for delta in np.linspace(-90,90,15):
        line = LineString(list(zip(lon+delta,lat)))
        ax.add_geometries([line],crs=geodetic,facecolor='none',edgecolor='k')

cartopy_example

Partial environment definition

Operating system

Mac OS X 10.15.4

Cartopy version

0.18.0

conda list

...
python                    3.7.6           h90870a6_5_cpython    conda-forge
cartopy                   0.18.0           py37h08e9697_0    conda-forge
shapely                   1.7.0            py37hfcf0db4_3    conda-forge
...
@QuLogic
Copy link
Member

QuLogic commented May 22, 2020

For Mollweide, I think this may be because of the fixed threshold:

cartopy/lib/cartopy/crs.py

Lines 1805 to 1807 in 60ae9a6

@property
def threshold(self):
return 1e5

Unfortunately, not very many projections scale the threshold with the globe.

QuLogic added a commit to QuLogic/cartopy that referenced this issue Jul 12, 2020
Fixes SciTools#1572, testing the example in the other projections as well.
QuLogic added a commit to QuLogic/cartopy that referenced this issue Jul 13, 2020
Fixes SciTools#1572, testing the example in the other projections as well.
QuLogic added a commit to QuLogic/cartopy that referenced this issue Jul 13, 2020
Using the example from SciTools#1572 shows that the lines are jagged and the
threshold is too high.
QuLogic added a commit to QuLogic/cartopy that referenced this issue Aug 6, 2020
Fixes SciTools#1572, testing the example in the other projections as well.
QuLogic added a commit to QuLogic/cartopy that referenced this issue Aug 6, 2020
Using the example from SciTools#1572 shows that the lines are jagged and the
threshold is too high.
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