TopoTijdReis kaartlagen in QGIS

Aanpassing 13 oktober 2020: meer service-lagen t/m 2019 zijn nu beschikbaar, scriptje en projectfile aangepast

Aanpassing 3 september 2020: de WMTS service werkt nu alleen nog onder httpS. Let dus op de code is aangepast (de service heeft een dag of wat de verkeerde url’s (http) in de capabilities gehad, maar per 3sept OK)

Aanpassing 19 juni 2020: voor de hele ongeduldige mensen (of bij wie het stukje code niet lukt): hier een gezipt QGIS topotijdreisqgisproject. Even uitpakken en het .qgz bestand openen als project.

Aanpassing 15 april 2018: code geupdate, werkt nu in zowel QGIS2 als QGIS3

Al weer even geleden presenteerde Het Kadaster samen met ESRI the website TopoTijdReis. Een website met prachtig kaartmateriaal: verdeelt over 75 kaartlagen de scans van topografische data van Het Kadaster van 1815 tot en met 2015!

Een website is mooi, maar dezelfde data kunnen gebruiken in bijvoorbeeld QGIS is nog leuker. Toen iemand mij benaderde of dat mogelijk was, heb ik geprobeerd om de lagen toe te voegen als ‘ArcGIS Mapserver layer’, wat mogelijk zou moeten zijn met QGIS 2.16 en hoger. Maar… ik kreeg het niet voor elkaar (en ernaar ‘googelende’ blijkbaar meer mensen niet)…
Maar… na wat gemopper van mij op Twitter dat ‘ze’ toch eens wat meer dit soort zaken m.b.v. OpenStandaarden zouden moeten ontsluiten, kreeg ik van het Kadaster het vriendelijke antwoord, dat de kaartlagen ook als WMTS beschikbaar zouden moeten zijn. Dat wil ik wel zien.
Eerst maar eens proberen op de WMTS link te klikken in een pagina als
https://tiles.arcgis.com/tiles/nSZVuSZjHpEZZbRo/arcgis/rest/services/Historische_tijdreis_1950/MapServer…
Mmm, geen tiles. Gelukkig werkt het getCapabilities document wel:
En op basis daarvan is het in QGIS eenvoudig om de service te bevragen. Voer de capabilities url (of de verkorte versie zonder het /1.0.0/WMTSCapabilities.xml deel) in in het ‘WMS/WMTS kaartlaag toevoegen’ dialoog

en de dialoog toont welke kaartlaag er op te halen is.

Prachtig! Werkt!

Maar het mooie aan de TopoTijdReis is natuurlijk dat je je de historie van een omgeving kunt bekijken door even door te tijd te reizen… We willen natuurlijk al die 75 kaartlagen ‘even’ laden in QGIS. Jammer dat elke laag in een aparte service is gedefineerd 🙁
Ik zie twee mogelijkheden:
1) we zorgen dat die 75 lagen komen te staan in het pdok.json lagenbestand van de PDOKServicePlugin, maar dan nog moet je 75 keer klikken om alle lagen te bekijken (deze optie stel ik dus even uit tot iemand hier heel vriendelijk om vraagt…)
2) we schrijven wat Python code welke we in de Python Console kunnen uitvoeren om die 75 lagen te laden.
DOEN:
– open de Python Console die je vind als menu item in het Plugins menu
– open een bestandje door op het ‘nieuw bestand’ knopje te klikken (4e knopje van links)
– kopieer daarin de volgende regels:

OPMERKING: onderstaande code zou zowel in recente QGIS 2.x als QGIS 3.x versies moeten werken:

# create a layer group
layer_group = QgsProject.instance().layerTreeRoot().insertGroup(0, 'TopoTijdReis')
# make the group mutually exlusive (one layer visible at a time)
layer_group.setIsMutuallyExclusive(True)
# copied from https://tiles.arcgis.com/tiles/nSZVuSZjHpEZZbRo/arcgis/rest/services/
years = ['1815', '1820', '1821', '1823_1829', '1830_1849', '1850', '1857', '1858', '1861', '1862', '1865', '1866', '1868', '1870', '1871', '1872', '1880', '1883', '1886', '1888', '1889', '1893', '1899', '1900', '1901', '1902', '1904', '1908', '1909', '1910', '1912', '1915', '1918', '1919', '1920', '1922', '1925', '1929', '1931', '1935', '1937', '1940', '1942', '1943', '1947', '1948', '1949', '1950', '1951', '1952', '1953', '1955', '1962', '1963', '1965', '1970', '1971', '1973', '1975', '1976', '1978', '1980', '1984', '1988', '1990', '1994', '1995', '1996', '1997', '1999', '2000', '2001', '2002', '2003', '2004', '2005', '2006', '2007', '2008', '2009', '2010', '2011', '2012', '2013', '2014', '2015',
'2016', '2017', '2018', '2019']
# for every year, create a layer
for year in years:
  layer = QgsRasterLayer("tileMatrixSet=default028mm&crs=EPSG:28992&layers=Historische_tijdreis_{0}&styles=default&format=image/jpgpng&url=https://tiles.arcgis.com/tiles/nSZVuSZjHpEZZbRo/arcgis/rest/services/Historische_tijdreis_{1}/MapServer/WMTS/".format(year, year), "TopoTijdReis {0}".format(year), "wms")
  # add it to the layer registry, without showing it yet
  # False, meaning not ready to add to legend
  if hasattr(qgis.core, 'QgsMapLayerRegistry'):
    # QGIS 2.x
    QgsMapLayerRegistry.instance().addMapLayer(layer, False)
  else:
    # QGIS 3.x
    QgsProject.instance().addMapLayer(layer, False)
  # add it on top of the group
  layer_group.insertLayer(0, layer)

Door nu op het blauwe ‘Play’-knopje te klikken wordt bovenstaande code uitgevoerd, en moet QGIS er ongeveer zo uit zien:

Zoals je ziet staan de 75 beschikbare kaartlagen (die ik trouwens van de ‘services’ url heb http://tiles.arcgis.com/tiles/nSZVuSZjHpEZZbRo/arcgis/rest/services/ gehaald) netjes in een aparte groep.
En het mooi van een layer-group is dat je ‘m ‘mutual exlusive’ kunt maken, d.w.z. er is slechts 1 kaartlaag tegelijk zichtbaar: bij het zichtbaar maken van een andere kaartlaag wordt een evt andere zichtbaar kaartlaag uit die groep verborgen.
Eigenlijk hebben we zo onze eigen TopoTijdReis in QGIS.
En mocht je nu via de python console gewoon maar 1 kaartlaag willen laden, gebruik dan:
iface.addRasterLayer(
"tileMatrixSet=default028mm&crs=EPSG:28992&layers=Historische_tijdreis_1820&styles=default&format=image/jpgpng&url=https://tiles.arcgis.com/tiles/nSZVuSZjHpEZZbRo/arcgis/rest/services/Historische_tijdreis_1820/MapServer/WMTS/", # uri
"TopoTijdReis 1820", # name for layer (as seen in QGIS)
"wms"
)
Een paar leuke screenshots nog ter afsluiting:
Zoals je ziet is de kaarttekst van 1815 nog franstalig. Wanneer namen we ook al weer afscheid van de Fransen? De daaropvolgende kaart is dan ook opgedragen “aan zijne majesteit den koning”

Verder is het goed te zien dat er hier en daar wel flink ‘geplakt’ is.
Letterlijk (ik denk door het oude Kadaster):

Figuurlijk (door de cache bouwer):

Een paar tips:
– Gebruik het ‘Tile Scale Panel/Paneel Schaal tegel’ om de juiste kaartschaal te bepalen: maak een kaartlaag aktief/zichtbaar en gebruik dan de schuif om ‘m op de juiste schaal te krijgen. (zie de grijze ‘schuif’ aan de rechterkant in een van de volgende screenshots)
– Zorg dat je in Extra/Opties/Rendering het vinkje ‘Gelijktijdig rendereren van kaartlagen …’ aan hebt. Voor Tiled layers scheelt dit aanzienlijk en ziet er VEEL beter uit!
– Gebruik de middelmuisknop om netjes met een factor 2 in of uit te zoomen (gebruik diezelfde middelmuisknop trouwens ook eens om eenvoudig te ‘pannen’ zonder de pan-tool aktief te hoeven maken
En om toch nog wat te mopperen te hebben:
– waarom heeft het kadaster/esri zich eigenlijk niet gehouden aan de Nederlands Tiling richtlijn? De kaartschalen komen niet overeen met die van de overige PDOK kaartlagen?
– waarom zijn juist de modernste kaartlagen niet scherp?
Maar verder: bedankt Kadaster en ESRI voor het beschikbaarstellen van dit prachtige kaartmateriaal!!