"""Show the same XYZ basemaps in Python with folium, leafmap and contextily.
Writes two HTML maps and one PNG figure (built from QGIS renders plus a
contextily render) to ./output.
"""
from pathlib import Path
import contextily as cx
import folium
import leafmap.foliumap as leafmap
import matplotlib.pyplot as plt
import matplotlib.image as mpimg
import xyzservices.providers as xyz
from pyproj import Transformer
OUT = Path(__file__).parent / "output"
OUT.mkdir(exist_ok=True)
OSM_URL = "https://tile.openstreetmap.org/{z}/{x}/{y}.png"
OSM_ATTR = '© OpenStreetMap contributors'
EOX_URL = "https://tiles.maps.eox.at/wmts/1.0.0/s2cloudless_3857/default/g/{z}/{y}/{x}.jpg"
EOX_ATTR = ('Sentinel-2 cloudless - https://s2maps.eu '
"by EOX IT Services GmbH (contains modified Copernicus Sentinel data 2016)")
CENTRE = [53.75, 0.01] # Withernsea, Holderness coast
# folium: OpenStreetMap is the default; extra layers need url + attribution
m = folium.Map(location=CENTRE, zoom_start=11, tiles=None)
folium.TileLayer(OSM_URL, attr=OSM_ATTR, name="OpenStreetMap", max_zoom=19).add_to(m)
folium.TileLayer(EOX_URL, attr=EOX_ATTR, name="Sentinel-2 cloudless 2016", max_zoom=14).add_to(m)
folium.TileLayer(xyz.Esri.WorldImagery.build_url(), attr=xyz.Esri.WorldImagery.html_attribution,
name="Esri World Imagery").add_to(m)
folium.LayerControl().add_to(m)
m.save(OUT / "folium_basemaps.html")
print("folium", folium.__version__, "-> folium_basemaps.html")
# leafmap: named basemaps or any XYZ URL
lm = leafmap.Map(center=CENTRE, zoom=11)
lm.add_tile_layer(url=EOX_URL, name="Sentinel-2 cloudless 2016", attribution=EOX_ATTR)
lm.to_html(str(OUT / "leafmap_basemaps.html"))
print("leafmap -> leafmap_basemaps.html")
# contextily: static basemap under a matplotlib axis (EPSG:3857 metres)
west, south, east, north = -0.10, 53.70, 0.12, 53.80
fig, axes = plt.subplots(1, 3, figsize=(15, 3.9), constrained_layout=True)
for ax, fname, title in [
(axes[0], "qgis_osm.png", "QGIS 3.40 XYZ: OpenStreetMap"),
(axes[1], "qgis_eox_2016.png", "QGIS 3.40 XYZ: Sentinel-2 cloudless 2016"),
]:
ax.imshow(mpimg.imread(OUT / fname))
ax.set_title(title, fontsize=11)
ax.set_axis_off()
axes[0].text(0.99, 0.01, "© OpenStreetMap contributors", transform=axes[0].transAxes,
ha="right", va="bottom", fontsize=8, bbox=dict(fc="white", ec="none", alpha=0.8))
axes[1].text(0.99, 0.01, "Sentinel-2 cloudless 2016 by EOX (CC BY 4.0)\nContains modified Copernicus Sentinel data 2016",
transform=axes[1].transAxes, ha="right", va="bottom", fontsize=7,
bbox=dict(fc="white", ec="none", alpha=0.8))
to_merc = Transformer.from_crs("EPSG:4326", "EPSG:3857", always_xy=True)
x0, y0 = to_merc.transform(west, south)
x1, y1 = to_merc.transform(east, north)
axes[2].set_xlim(x0, x1)
axes[2].set_ylim(y0, y1)
cx.add_basemap(axes[2], crs="EPSG:3857", source=xyz.OpenStreetMap.Mapnik, zoom=11,
attribution=False,
headers={"User-Agent": "TechnicalGIS-tutorial/1.0 (+https://technicalgis.com)"})
print("contextily", cx.__version__, "extent (EPSG:3857 m)", [round(v) for v in (x0, y0, x1, y1)])
axes[2].set_title("contextily: OpenStreetMap.Mapnik, zoom 11", fontsize=11)
axes[2].set_axis_off()
axes[2].text(0.99, 0.01, "© OpenStreetMap contributors", transform=axes[2].transAxes,
ha="right", va="bottom", fontsize=8, bbox=dict(fc="white", ec="none", alpha=0.8))
fig.savefig(OUT / "basemaps_qgis_and_python.png", dpi=110)
print("saved basemaps_qgis_and_python.png")