Coordinates_operations — CRS transformations

This module provides coordinate transformation utilities between Belgian CRS:

Lambert 72 (EPSG:31370) — legacy Belgian grid
Lambert 2008 (EPSG:3812) — current Belgian grid
WGS 84 (EPSG:4326) — lat/lon

It also supports raster reprojection via GDAL.

[1]:

from wolfhece.Coordinates_operations import transform_coordinates
import numpy as np

Transforming point coordinates

transform_coordinates converts an (n, 2) numpy array from one EPSG to another. Large arrays are processed in chunks with multiprocessing.

[2]:

# Points in Lambert 2008 (EPSG:3812)
pts_3812 = np.array([
    [710000.0, 600000.0],  # approximate Liège area
    [650000.0, 650000.0],  # approximate Brussels area
])

# Convert to Lambert 72 (EPSG:31370)
pts_31370 = transform_coordinates(
    pts_3812,
    inputEPSG='EPSG:3812',
    outputEPSG='EPSG:31370'
)

print("Lambert 2008 -> Lambert 72:")
for a, b in zip(pts_3812, pts_31370):
    print(f"  ({a[0]:.1f}, {a[1]:.1f}) -> ({b[0]:.1f}, {b[1]:.1f})")

Lambert 2008 -> Lambert 72:
  (710000.0, 600000.0) -> (209992.5, 99991.9)
  (650000.0, 650000.0) -> (149997.9, 149999.4)

[3]:

# Convert to WGS 84 (lat/lon)
pts_wgs84 = transform_coordinates(
    pts_3812,
    inputEPSG='EPSG:3812',
    outputEPSG='EPSG:4326'
)

print("Lambert 2008 -> WGS 84 (lon, lat):")
for a, b in zip(pts_3812, pts_wgs84):
    print(f"  ({a[0]:.1f}, {a[1]:.1f}) -> ({b[0]:.6f}, {b[1]:.6f})")

Lambert 2008 -> WGS 84 (lon, lat):
  (710000.0, 600000.0) -> (5.209195, 50.207972)
  (650000.0, 650000.0) -> (4.368720, 50.660611)

Raster reprojection

reproject_and_resample_raster reprojects a GeoTIFF between CRS using GDAL.

from wolfhece.Coordinates_operations import reproject_and_resample_raster
from osgeo import gdal

reproject_and_resample_raster(
    input_raster_path='input_L2008.tif',
    output_raster_path='output_L72.tif',
    input_srs='EPSG:3812',
    output_srs='EPSG:31370',
    resampling_method=gdal.GRA_Bilinear,
    xRes=1.0, yRes=1.0
)

Summary

Function	Use case
`transform_coordinates(pts, inputEPSG, outputEPSG)`	Convert point arrays between CRS
`reproject_and_resample_raster(...)`	Reproject GeoTIFF rasters via GDAL

Common Belgian EPSG codes:

EPSG	Name
31370	Belgian Lambert 72
3812	Belgian Lambert 2008
4326	WGS 84 (lat/lon)