Hoe vraag je de hoogte op d.m.v. lat-long coordinaten?

Datasets AHN
,
#1

Gegeven een .TIF bestand van een raster verkregen via PDOK - AHN3 downloads, hoe kan je met behulp van een Python script de hoogte bepalen door middel van coordinaten?

Wat ik heb geprobeerd is hetvolgende:

from PIL import Image
from PIL.TiffTags import TAGS

with Image.open(r"M5_44CN2.TIF") as img:
    meta_dict= {TAGS[key] : img.tag[key] for key in img.tag.keys()}

meta_dict bevat een aantal keys maar geen van deze keys lijkt ook maar enige informatie te geven over de coördinaten van alle punten. Wel is het waar dat wanneer ik een ander raster kies, de key ‘ModalTiepointTag’ een andere x en y coordinaat bevat. Echter, is mij de referentie van deze coördinaten volstrekt onduidelijk.

Kan iemand uitleg geven over de informatie die het .tif bestand bevat? Specifiek, hoe kan ik de exacte locatie van de hoogtepunten uitlezen?

EDIT: Het lijkt erop dat de x en y in ‘ModelTiepointTag’ zogeheten Rijksdriehoek coördinaten zijn. Is dit juist, en slaat het dan op het centrum van het raster of een van de hoekpunten?

#2

In de header van de GeoTIFF staan de begrenzingen van het grid vermeld. Als ik dat met het commando gdalinfo M5_44CN2.TIF uitlees krijg ik:

Driver: GTiff/GeoTIFF
Files: M5_44CN2.TIF
Size is 1000, 1250
Coordinate System is:
PROJCRS["Amersfoort / RD New",
    BASEGEOGCRS["Amersfoort",
        DATUM["Amersfoort",
            ELLIPSOID["Bessel 1841",6377397.155,299.1528128,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4289]],
    CONVERSION["RD New",
        METHOD["Oblique Stereographic",
            ID["EPSG",9809]],
        PARAMETER["Latitude of natural origin",52.1561605555556,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",5.38763888888889,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9999079,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",155000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",463000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["easting (X)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["northing (Y)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Engineering survey, topographic mapping."],
        AREA["Netherlands - onshore, including Waddenzee, Dutch Wadden Islands and 12-mile offshore coastal zone."],
        BBOX[50.75,3.2,53.7,7.22]],
    ID["EPSG",28992]]
Data axis to CRS axis mapping: 1,2
Origin = (105000.000000000000000,412500.000000000000000)
Pixel Size = (5.000000000000000,-5.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  105000.000,  412500.000) (  4d39'51.40"E, 51d41'59.89"N)
Lower Left  (  105000.000,  406250.000) (  4d39'54.63"E, 51d38'37.63"N)
Upper Right (  110000.000,  412500.000) (  4d44'11.79"E, 51d42' 1.41"N)
Lower Right (  110000.000,  406250.000) (  4d44'14.70"E, 51d38'39.16"N)
Center      (  107500.000,  409375.000) (  4d42' 3.13"E, 51d40'19.54"N)
Band 1 Block=256x256 Type=Float32, ColorInterp=Gray
  NoData Value=3.4028234663852886e+38

Op basis van de Upper/Lower Left/Right kan je ieder pixel vinden. De linkerbovenhoek van de eerste pixel heeft RD-coördinaten x=105000 m; y=412500 m, de linkerbovenhoek van de tweede pixel heeft RD-coördinaten x=105005 m; y=412500 m, enz.

#3

PS: Als je niet tussen RD-coördinaten en lat-long wil omrekenen zijn er ook een WFS- en WCS-service waarmee je het AHN direct voor een bbox (van eventueel 1 pixel) in lat-long kunt bevragen: https://www.pdok.nl/…/actueel-hoogtebestand-nederland…

#4

Mogelijk kan rasterio dit werk ook wat vereenvoudigen (specifiek .sample()).