Luchtfoto's ophalen in R obv coördinaten

Datasets Luchtfoto's
#1

Hallo,

Mijn naam is Esmee Kramer en ik werk als data scientist bij de Omgevingsdienst Midden-West Brabant. Samen met een collega van de Omgevingsdienst Brabant Noord, ben ik een image recognition project gestart om natte koeltorens op te sporen. Wij hebben weinig kennis van GIS of andere spatial data, maar werken veel met R.

Als startpunt willen we de luchtfoto’s ophalen van de koeltorens die bij ons bekend zijn ophalen. We denken dat we het beste met de WMTS service kunnen werken. We hebben hier de adressen en coördinaten van (en kunnen we zelf omzetten naar een ander coördinatenstelsel). We proberen middels het tweaken van de volgende link de juiste luchtfoto’s naar voren te krijgen:

https://geodata.nationaalgeoregister.nl/luchtfoto/infrarood/wmts/2016_ortho25IR/EPSG:28992/{TileMatrix}/{TileCol}/{TileRow}.jpeg

Inmiddels hebben we gevonden dat TileMatrix voor het zoomniveau staat, TileCol voor (een equivalent van) het X-coördinaat, en TileRow voor (een equivalent van) het Y-coördinaat. Het lukt ons echter niet om de juiste luchtfoto’s te vinden. Een voorbeeld locatie is; 51.7008, 4.9416 (lat/lon). Hoe weten we welke waarden we bij TileCol en TileRow moeten invullen, en komen we tot de juiste url? Kan iemand ons daarbij helpen?

Bij voorbaat dank!
Met vriendelijke groet, Esmee

#2

https://wiki.openstreetmap.org/wiki/Slippy_map_tilenames

1 like
#3

Bedankt voor je reactie. Als ik hier de coördinaten 51.7008 en 4.9416 invul met bijv. zoom niveau 3, komt daaruit: 4, 2
Als ik deze invul in de link kom ik uit op:
https://geodata.nationaalgeoregister.nl/luchtfoto/infrarood/wmts/2016_ortho25IR/EPSG:4326/3/4/2.jpeg
De gezochte locatie (Waspik) ligt daar echter niet in … Kun je ons verder helpen?

#4

Hoi @ekrameromwb
De coordinaten die jij hebt zijn in graden mogelijk CRS:84 (EPSG:4326…?), in andere woorden in een bepaalde projectie.

De URL (dit jij wilt gebruiken?)

‘Spreekt’ op dit end-point alleen EPSG:28992, je zou uit de WMTSCapabilities het endpoint kunnen halen wat wel in om kan gaan met een projectie die graden spreekt (i.p.v. meters)

of de punten die je hebt omzetten naar EPSG:28992 met iets als gdaltransform

gdaltransform -s_srs EPSG:4326 -t_srs EPSG:28992
4.9416 51.7008
124195.740542834 412542.303586653 0
#5

Hoi Wouter,

Ik wil inderdaad graag werken met URL(s). Ik kan inderdaad EPSG:28992 aanpassen naar 4326 (lat/long) of 3857(met meter als eenheid). De luchtfoto in mijn vorige post met Flevoland er deels op, is gebaseerd op EPSG:4326, maar werkt bij mij dus blijkbaar niet helemaal naar behoren. Weet jij wat ik verkeerd doe?

#6

Hoi @ekrameromwb

Op zoomlevel 12, komt er bij mij 2104, 1358 uit
Als ik de url aanpas naar EPSG:3857, dan zit Waspik er wel bij.

https://geodata.nationaalgeoregister.nl/luchtfoto/rgb/wmts/2016_ortho25IR/EPSG:3857/12/2104/1358.jpeg

Dat is ook de url die ik zelf gebruik voor satellietbeelden van nederland (zonder de IR).
Als je 2016 veranderd naar 2018, dan heb je meer recente foto’s

Groet,
Hans

1 like
#7

Nee, verder dat je coordinaten in EPSG:3857 hebt, iets met luchtofotos wil doen en R gebruikt weet ik niet goed waar in die ‘combinatie’ mogelijk een fout zit. :smiley: Maar ik denk dat je al wel redelijk op de goede weg bent/zit.

Ik kan nu niet anders zeggen hoe ik dit issue zou (mogelijk) tackelen met de eerste stappen:

  1. keuze maken in welke projectie je zou willen werken, dit gezien mogelijk issue met herprojectie (lat/lon → lon/lat, meters, source van de brondata,… enz…)
  2. Mogelijk aan de hand van de ‘gegevens’ in de capabitilies: https://geodata.nationaalgeoregister.nl/luchtfoto/infrarood/wmts/1.0.0/WMTSCapabilities.xml de ‘juiste’ tegel berekenen. Door middel van initieel een sanity check

We weten dat je nu het punt hebt:

dat ligt ongeveer hier,

image
image927×684 214 KB

Dus dat zou betekenen dat het WMTS tegel …/EPSG:3857/17/67335/43472.jpeg zou betreffen

image

Mogelijk dat je in de R code die je hebt iets kan delen, qua logica of packages die je gebruikt?
Dat je met de ‘uiteindelijke’ tegel die je zou moeten hebben ‘terug kan redeneren’?

1 like
#8

Als ik de gegeven coordinaten invoer in Google maps, dan komt daar Industrieweg 3 in Waspik uit, bij het bedrijf aldaar, Maasoever Agriprocessing, zie ik koelers op het dak staan. Ik denk dat Esmee dat adres bedoelt.
Dus ergens zit er een verschuiving van coordinaten, een fout in de berekening, of een conversie fout van de ene projectie naar de andere.
Ik heb geen ervaring met R, maar ik ga eens kijken of ik de juiste tile kan vinden. Als je weet naar welke url je toe moet werken, zie je de fout misschien.

1 like
#9

Ik heb de antwoorden hierboven wat samengevoegd en toen het een en ander geprobeerd. Ik ben ervan uitgegaan dat het adres dat ik hier boven genoemd heb de locatie is die je zoekt.

Let op het jaar, in de url die je gepost hebt, gebruik je de data uit 2016. Er is ook data van 2018 en 2019. Mooi is om te zien dat er een andere kast op het dak van het bovengenoemde adres staat in 2019 dan 2018. Je kan het jaar ook vervangen door Actueel (let op de hoofdletter) om altijd de meest recente foto te hebben. Ik heb wel gemerkt dat er wel kwaliteits verschil zit tussen de verschillende jaren.

Ik heb eerst de voorbeeldcode gepakt van de slippy map tilenames beschrijving uit de link van leptonix hierboven. Ikheb c# gebruikt omdat de omgevng is waar ik het meest in werk. Daarna de projectie aangepast, zoals leptonix aangeeft. Dan kom ik op de volgende urls

rgb:
zoom 12 => https://geodata.nationaalgeoregister.nl/luchtfoto/rgb/wmts/Actueel_ortho25/EPSG:3857/12/2104/1358.jpeg

zoom 19 => https://geodata.nationaalgeoregister.nl/luchtfoto/rgb/wmts/Actueel_ortho25/EPSG:3857/19/269340/173885.jpeg

ir:

zoom 12 => https://geodata.nationaalgeoregister.nl/luchtfoto/infrarood/wmts/Actueel_ortho25IR/EPSG:3857/12/2104/1358.jpeg

zoom 19 => https://geodata.nationaalgeoregister.nl/luchtfoto/infrarood/wmts/Actueel_ortho25IR/EPSG:3857/19/269340/173885.jpeg

Als je voorbeelden van andere zoom levels wilt hebben, laat het dan even weten. Graag ook even laten weten of het gelukt is.

1 like
#10

Bedankt voor jullie reacties en voor het meedenken, het is nu gelukt!

Er zat inderdaad een mis match in de projecties. We gebruiken nu de oorspronkelijke latitude & longitude waarden, en daar hoort dan de volgende url bij (zoals jij ook beschrijft @EwoutdeBoer, bedankt daarvoor!) :

“https://geodata.nationaalgeoregister.nl/luchtfoto/rgb/wmts/2018_ortho25/EPSG:3857/19/269342/173885.jpeg”

(19 is volgens mij het maximale zoomniveau)

De laatste 2 waarden in deze url kunnen we ophalen via de volgende functie die op wiki pagina stond (thanks @leptonix);

deg2num <-function(lat_deg, lon_deg, zoom){
lat_rad ← lat_deg * pi /180
n ← 2.0 ^ zoom
xtile ← floor((lon_deg + 180.0) / 360.0 * n)
ytile = floor((1.0 - log(tan(lat_rad) + (1 / cos(lat_rad))) / pi) / 2.0 * n)
return( c(xtile, ytile))
}

Inmiddels hebben we dit geautomatiseerd kunnen doen voor gelijk alle 8 omliggende luchtfoto’s (zodat je een raster van 3x3 foto’s ophaalt, voor het geval een koeltoren niet in het geheel om de eerste, middelste foto staat).

Nu op naar de volgende stap van labeling!

2 likes
#11