check scaling of basemaps and axis

[veenstrajelmer]

I noticed that the contextily basemap scales my figure axes significantly more than a geopandas object does. This is documented for my tool here: Deltares/dfm_tools#966. Is this something that you recognize or is this even intented? I am curious to hear. If it is unexpected I can share an example with less dependencies.

[martinfleis]

Would you be able to provide a reproducible example? GeoPandas sets aspect to equal for projected CRS and to custom for geographical. Contextily does not change the aspect afaik but it does assume equal due to pixels of the underlying tile image being square. Could be a result of some interplay between these but a minimal example would help here.

[veenstrajelmer]

Yes, I created one:

import matplotlib.pyplot as plt
plt.close('all')
import contextily as ctx
from contextily.plotting import _reproj_bb
import geopandas as gpd

def plot_coastlines(ax, crs=None):
    xlim = ax.get_xlim()
    ylim = ax.get_ylim()
    xmin, xmax = xlim
    ymin, ymax = ylim
    # xmin, xmax, ymin, ymax = -180, 180, -90, 90
    if crs is not None:
        xmin, xmax, ymin, ymax = _reproj_bb(xmin, xmax, ymin, ymax, 
                                            s_crs=crs, t_crs='EPSG:4326')

    # gshhg-shp-2.3.7.zip can be retrieved from https://www.ngdc.noaa.gov/mgg/shorelines/data/gshhg/latest/
    file_shp = r"c:\Users\veenstra\Downloads\gshhg-shp-2.3.7\GSHHS_shp\c\GSHHS_c_L1.shp"
    bbox = (xmin, ymin, xmax, ymax)
    coastlines_gdb = gpd.read_file(file_shp, columns=[], bbox=bbox)
    if crs:
        coastlines_gdb = coastlines_gdb.to_crs(crs)
    
    coastlines_gdb.plot(ax=ax, edgecolor='k', facecolor='none', zorder=5)
    ax.set_xlim(xlim)
    ax.set_ylim(ylim)

lon_min, lon_max, lat_min, lat_max = -6, 3.5, 48, 51.5 # france

plot_crs = 'EPSG:4326'
# plot_crs = 'EPSG:3857'

lon_min, lon_max, lat_min, lat_max = _reproj_bb(lon_min, lon_max, lat_min, lat_max, s_crs='EPSG:4326', t_crs=plot_crs)

fig, (ax1, ax2) = plt.subplots(2,1)
ax1.set_xlim(lon_min, lon_max)
ax1.set_ylim(lat_min, lat_max)
ctx.add_basemap(ax=ax1, crs=plot_crs, attribution=False)
plot_coastlines(ax=ax1, crs=plot_crs)

ax2.set_xlim(lon_min, lon_max)
ax2.set_ylim(lat_min, lat_max)
plot_coastlines(ax=ax2, crs=plot_crs)
ctx.add_basemap(ax=ax2, crs=plot_crs, attribution=False)

With plot_crs = 'EPSG:4326' (different aspect ratio):

With plot_crs = 'EPSG:3857' (equal aspect ratio):

With plot_crs = 'EPSG:4326' and xmin, xmax, ymin, ymax = -180, 180, -90, 90 (very comparable aspect ratio):

It seems that the scaling is different depending on the bounding box of the object to be plotted. When I set the full globe as bounding box for the coastline plotting, the axis scaling gets close to the one by contextily. I therefore assumed that the bounding box for contextily images is always the full globe, but I do see that the axis limits are taken into account in the warping procedure:

contextily/contextily/plotting.py

Lines 112 to 139 in 0e38541

	xmin, xmax, ymin, ymax = ax.axis()
	if isinstance(source, str):
	try:
	source = providers.query_name(source)
	except ValueError:
	pass
	# If web source
	if (
	source is None
	or isinstance(source, (dict, TileProvider))
	or (isinstance(source, str) and source[:4] == "http")
	):
	# Extent
	left, right, bottom, top = xmin, xmax, ymin, ymax
	# Convert extent from `crs` into WM for tile query
	if crs is not None:
	left, right, bottom, top = _reproj_bb(
	left, right, bottom, top, crs, "epsg:3857"
	)
	# Download image
	image, extent = bounds2img(
	left, bottom, right, top, zoom=zoom, source=source, ll=False, zoom_adjust=zoom_adjust
	)
	# Warping
	if crs is not None:
	image, extent = warp_tiles(image, extent, t_crs=crs, resampling=resampling)

So I do not understand yet why the behaviour is different. I am probably missing something, so please feel free to clarify.

The shapefile with coastlines:
gshhg-shp-2.3.7_coarse_only.zip

[martinfleis]

Thanks! What happens here is that when using EPSG 4326, GeoPandas tried so set the aspect to minimise distortion of geometries which are represented in degrees on a sphere. And contextily resets the aspect to 1:1. We should probably avoid doing that but given it is not done explicitly, not sure how easy that is.

[veenstrajelmer]

Amazing, thanks for fixing this!