[GRASS-dev] shift raster when projecting

Dear all,

I am trying to project a DSM raster layer from Belgium (https://overheid.vlaanderen.be/en/producten-diensten/digital-elevation-model; EPSG 31370 to EPSG 28992 (Amersfoort / RD New).

The result is a DEM raster layer with the features shifted roughly 100 m as compared to the same features in the AHN (Dutch DEM, https://www.ahn.nl/ahn-4). See screenshot and code below.

Reprojecting the same in QGIS using the warp function (gdalwarp) results in the Belgium DSM layer to be fairly closely aligned to the Dutch one.

I tried to set the bounds before, align to the target raster, and other suggestions in the help file of r.proj, to no avail. Any idea what I might be doing wrong?

Cheers,

Paulo

image795.png

AHN (Netherlands)

EPSG 28992
projection: 99 (Amersfoort / RD New)
zone: 0
datum: rd18
ellipsoid: bessel
north: 397562.15979713
south: 373497.15979713
west: 107708.49748491
east: 133913.49748491
nsres: 5
ewres: 5
rows: 4813
cols: 5241
cells: 25224933

Belgium DSM

EPSG 31370
projection: 99 (BD72 / Belgian Lambert 72)
zone: 0
datum: bel72
ellipsoid: international
north: 250000
south: 238000
west: 162000
east: 194000
nsres: 5
ewres: 5
rows: 2400
cols: 6400
cells: 15360000

Reproject

(Thu Jun 15 17:11:08 2023)
r.proj location=Merkske_Belgium mapset=PERMANENT input=DHMVIIDSMRAS5m_k02 output=test memory=3000 resolution=5
Selected PROJ pipeline:
+proj=pipeline +step +inv +proj=sterea +lat_0=52.1561605555556 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +step +proj=lcc +lat_0=90 +lon_0=4.36748666666667 +lat_1=51.1666672333333 +lat_2=49.8333339 +x_0=150000.013 +y_0=5400088.438 +ellps=intl


Selected PROJ pipeline:
+proj=pipeline +step +inv +proj=lcc +lat_0=90 +lon_0=4.36748666666667 +lat_1=51.1666672333333 +lat_2=49.8333339 +x_0=150000.013 +y_0=5400088.438 +ellps=intl +step +proj=sterea +lat_0=52.1561605555556 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel


Input:
Cols: 4114 (original: 6400)
Rows: 2400 (original: 2400)
North: 250000.000000 (original: 250000.000000)
South: 238000.000000 (original: 238000.000000)
West: 173430.000000 (original: 162000.000000)
East: 194000.000000 (original: 194000.000000)
EW-res: 5.000000
NS-res: 5.000000

Output:
Cols: 4108 (original: 5241)
Rows: 2456 (original: 4813)
North: 396852.159797 (original: 397562.159797)
South: 384572.159797 (original: 373497.159797)
West: 107708.497485 (original: 107708.497485)
East: 128248.497485 (original: 133913.497485)
EW-res: 5.000000
NS-res: 5.000000

Allocating memory and reading input raster map…
Selected PROJ pipeline:
+proj=pipeline +step +inv +proj=sterea +lat_0=52.1561605555556 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +step +proj=lcc +lat_0=90 +lon_0=4.36748666666667 +lat_1=51.1666672333333 +lat_2=49.8333339 +x_0=150000.013 +y_0=5400088.438 +ellps=intl


On Fri, Jun 16, 2023 at 7:39 AM Paulo van Breugel
<p.vanbreugel@gmail.com> wrote:

Dear all,

I am trying to project a DSM raster layer from Belgium (https://overheid.vlaanderen.be/en/producten-diensten/digital-elevation-model; EPSG 31370 to EPSG 28992 (Amersfoort / RD New).

The result is a DEM raster layer with the features shifted roughly 100 m as compared to the same features in the AHN (Dutch DEM, https://www.ahn.nl/ahn-4). See screenshot and code below.

Reprojecting the same in QGIS using the warp function (gdalwarp) results in the Belgium DSM layer to be fairly closely aligned to the Dutch one.

I tried to set the bounds before, align to the target raster, and other suggestions in the help file of r.proj, to no avail. Any idea what I might be doing wrong?

Just guessing, maybe some geodetic datum file missing?

Maybe try this?

# https://proj.org/usage/environmentvars.html#envvar-PROJ_NETWORK
export PROJ_NETWORK=ON

Markus