cartopy tiles in Orthographic projection miss some datas

[ludwigVonKoopa]

Description

When using tiles in Orthographic projection, i noticed missing data.

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from cartopy.io.img_tiles import Stamen

fig = plt.figure(figsize=(5,5), dpi=120)
ax = fig.add_subplot(1, 1, 1, projection=ccrs.Orthographic(-10, 60))
tiler = Stamen('terrain-background')

ax.add_image(tiler, 3)

I noticed white data at the top of the globe. It doesn't appear with precision < 3

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from cartopy.io.img_tiles import Stamen

tiler = Stamen('terrain-background')
fig = plt.figure(figsize=(5,5), dpi=120)
for x in range(1,5):
    # print(x)
    ax = fig.add_subplot(2, 2, x, projection=ccrs.Orthographic(-10, 60))

    ax.add_image(tiler, x)
    ax.set_title(f"tile pre={x}")

Is that something wanted ?

Thanks

Cartopy version

import cartopy
cartopy.__version__

0.20.1

[greglucas]

It looks like the stamen-terrain doesn't go fully up to the pole. There is a band near the pole in PlateCarree as well.

ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())
ax.set_global()

I think this is unactionable on Cartopy's side.

[QuLogic]

I don't think this is referring to the pole, but the top of the circle; one tile that only shows a few pixels is not added.

[ludwigVonKoopa]

I checked this issue, i found the "problem".

When cartopy load the image, it load only the tiles who will be seen in the final images (i.e. skip tiles which are behind the earth).
This is the PlateCarree representation of the projected tiles loaded, it is more understandable :

code to create figure

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from cartopy.io.img_tiles import Stamen
tiler = Stamen(style="terrain-background", cache="/tmp/tiles/")


fig = plt.figure(figsize=(10,8), dpi=120)

proj = ccrs.Orthographic(0, 45)

ax_pla = fig.add_subplot(2, 1, 1, projection=ccrs.PlateCarree(), frameon=True)
ax_ort = fig.add_axes([0., 0.32, 1, 0.3], projection=proj)

# this is what cartopy.mpl.geoaxes.GeoAxes.draw do for tiles computing
geom = ax_ort._get_extent_geom(tiler.crs)
img, extent, origin = tiler.image_for_domain(geom, 4)

for ax in [ax_pla, ax_ort]:
    ax.imshow(img, extent=extent, origin=origin, transform=tiler.crs)
    
ax_pla.set_title("tiles loaded for ccrs.Orthographic(0, 45) projection")

We see why there is no water sometimes on the border of the earth : tiles are not loaded.

To have a better representation, i plot the projection border on the plot :

code to create figure

import matplotlib.patches as mpatches
import shapely.geometry as sgeom

bounds_proj = sgeom.Polygon(proj.boundary)

for ax in [ax_pla, ax_ort]:
    c_bounds = ax.add_geometries([bounds_proj], crs=proj, color="green", alpha=0.3, zorder=40)
    
geo2 = mpatches.Patch(color='green', label="ortho boundaries")
ax_pla.legend(handles=[geo2], loc='lower right')

We see white patch inside the green polygon => it should not.
But why are those tiles not loaded ?

Because the boundaries computed are different from the theorical boundaries :

code to create figure

for ax in [ax_pla, ax_ort]:
    c_geom = ax.add_geometries([geom], crs=tiler.crs, color="red", alpha=0.3, zorder=30)
    
geo1 = mpatches.Patch(color='red', label='eroded boundaries')
ax_pla.legend(handles=[geo1, geo2], loc='lower right')

Red boundaries correspond to tiles loaded.
This eroded boundaries are computed in the _get_extent_geom function with the .threshold attribute.

If we use a lower threshold, we see that there is no more problems

fig = plt.figure(figsize=(10,8), dpi=120)

proj = ccrs.Orthographic(0, 45)
proj.threshold /= 10  # <= lower threshold

ax_pla = fig.add_subplot(2, 1, 1, projection=ccrs.PlateCarree(), frameon=True)
ax_ort = fig.add_axes([0., 0.32, 1, 0.3], projection=proj)

bounds_proj = sgeom.Polygon(proj.boundary)
geom = ax_ort._get_extent_geom(tiler.crs)
img, extent, origin = tiler.image_for_domain(geom, 4)

for ax in [ax_pla, ax_ort]:
    ax.imshow(img, extent=extent, origin=origin, transform=tiler.crs)
    c_bounds = ax.add_geometries([bounds_proj], crs=proj, color="green", alpha=0.3, zorder=40)
    c_geom = ax.add_geometries([geom], crs=tiler.crs, color="red", alpha=0.3, zorder=30)
    
ax_pla.set_title("tiles loaded for ccrs.Orthographic(0, 45) projection, ")

geo2 = mpatches.Patch(color='green', label="ortho boundaries")
geo1 = mpatches.Patch(color='red', label='eroded boundaries')
ax_pla.legend(handles=[geo1, geo2], loc='lower right')

In the code, it is stated that the boundaries are eroded "to avoid transform issues". What are the cases of those problems ?

So i ended with this conclusion :

• The threshold value is too big, but
• This threshold problem may be too specific, and don't justify a threshold modification (mainly because it can be solved with proj.threshold /= 10 and
• threshold modification may introduce transform issues

What do you think about this ? Should we modify the threshold value ?

Thanks

[greglucas]

Great investigation! Nice figures to show exactly what is going on.

It looks like the boundary calculation is already pre-eroded in the CRS as well and that could potentially be re-addressed too.
#1154 (comment)

I agree with your threshold issue though, right now it seems like it is too large. It is about a factor of 10 larger (relatively) than the PlateCarree threshold of 0.5. So, it seems like perhaps rather than the 2% factor, we could reduce that to something closer to 0.001.

cartopy/lib/cartopy/crs.py

Line 2074 in 0f60502

self.threshold = np.diff(self._x_limits)[0] * 0.02

I think the boundary erosion issue would likely take a lot of looking at corner cases and figuring out why tests are failing. It would be nice to remove that, but I think that would be a pretty big undertaking, so I would focus on the threshold to fix this issue.

[ludwigVonKoopa]

Ok i agree too to fix the threshold, i'll submit a PR.

p.s.: how did you put the code line this way ? it quite pretty!

Thanks

[greglucas]

👍
If you click on the line number within the code (or hold shift to get multiple lines), there is a "copy permalink" button and then just paste that link in the issue and GitHub will do the nice rendering.

[ludwigVonKoopa]

Oh! Thanks you!