Calculate the Depths of Offshore Wind Turbines in the North Sea — GDAL-only Reproduction¶

This notebook reproduces the workflow of _1_determine_offshore_turbine_depths.ipynb using the GDAL/OGR Python bindings directly to compare it geokit.

1. Data Preparation and Script Setup¶

1.1 Imports¶

We use osgeo.gdal and osgeo.ogr for all geospatial operations, plus numpy for array arithmetic and matplotlib for plotting. ZenodoDataDownloader is reused only as a download helper — it is not part of the geoprocessing pipeline.

In [ ]:

loc-setup

import pathlib

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Polygon as MplPolygon
from matplotlib.collections import PatchCollection
from osgeo import gdal, ogr, osr

from geokit.get_test_data import ZenodoDataDownloader  # download helper only

gdal.UseExceptions()
ogr.UseExceptions()

current_working_directory = pathlib.Path.cwd()

import pathlib import numpy as np import matplotlib.pyplot as plt from matplotlib.patches import Polygon as MplPolygon from matplotlib.collections import PatchCollection from osgeo import gdal, ogr, osr from geokit.get_test_data import ZenodoDataDownloader # download helper only gdal.UseExceptions() ogr.UseExceptions() current_working_directory = pathlib.Path.cwd()

1.2 Download the input data¶

Same three datasets as in the GeoKit example: IHO sea basins (shapefile), 2050 offshore wind turbine locations (shapefile), and pre-tiled GEBCO bathymetry (GeoTIFFs).

In [ ]:

loc-download

iho_download = (
    "https://geo.vliz.be/geoserver/wfs?request=getfeature&service=wfs&version=1.0.0&typename=MarineRegions:iho&outputformat=SHAPE-ZIP",
    "iho.zip",
    "iho",
    {},
)
turbine_location_download = (
    "https://zenodo.org/record/10259046/files/turbine_locations shapefile.zip?download=1",
    "turbine_locations_shapefile.zip",
    "turbine_locations_shapefile",
    {},
)
gebco_tiles_download = (
    "https://zenodo.org/record/17047388/files/gebco_tiles.zip?download=1",
    "gebco_tiles.zip",
    "gebco_tiles",
    {},
)

data_downloader = ZenodoDataDownloader()
(
    path_to_sea_basins_folder,
    path_to_offshore_wind_points_folder,
    path_to_bathymetry_rasters_folder,
) = data_downloader.download_and_extract_parallel(
    download_list=[iho_download, turbine_location_download, gebco_tiles_download]
)

path_to_sea_basins_shape_file = str(path_to_sea_basins_folder.joinpath("iho.shp"))
offshore_wind_points_path_str = str(path_to_offshore_wind_points_folder)
bathymetry_rasters_dir_str = str(path_to_bathymetry_rasters_folder)

iho_download = ( "https://geo.vliz.be/geoserver/wfs?request=getfeature&service=wfs&version=1.0.0&typename=MarineRegions:iho&outputformat=SHAPE-ZIP", "iho.zip", "iho", {}, ) turbine_location_download = ( "https://zenodo.org/record/10259046/files/turbine_locations shapefile.zip?download=1", "turbine_locations_shapefile.zip", "turbine_locations_shapefile", {}, ) gebco_tiles_download = ( "https://zenodo.org/record/17047388/files/gebco_tiles.zip?download=1", "gebco_tiles.zip", "gebco_tiles", {}, ) data_downloader = ZenodoDataDownloader() ( path_to_sea_basins_folder, path_to_offshore_wind_points_folder, path_to_bathymetry_rasters_folder, ) = data_downloader.download_and_extract_parallel( download_list=[iho_download, turbine_location_download, gebco_tiles_download] ) path_to_sea_basins_shape_file = str(path_to_sea_basins_folder.joinpath("iho.shp")) offshore_wind_points_path_str = str(path_to_offshore_wind_points_folder) bathymetry_rasters_dir_str = str(path_to_bathymetry_rasters_folder)

2. Load, Extract and Display the Region of the North Sea¶

GeoKit uses geokit.vector.extractFeatures(..., where="NAME='North Sea'") which under the hood calls OGRLayer::SetAttributeFilter. We do the same here directly through OGR. The geometry is cloned so it survives after the data source is closed.

2.1 Load and Extract the Region of the North Sea¶

In [3]:

loc-workflow

EPSG4326 = osr.SpatialReference()
EPSG4326.ImportFromEPSG(4326)
# Force lon/lat axis ordering for consistency with GDAL <3.x and the GeoKit example.
EPSG4326.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)

iho_ds = ogr.Open(path_to_sea_basins_shape_file)
iho_layer = iho_ds.GetLayer()
iho_layer.SetAttributeFilter("NAME = 'North Sea'")

north_sea_feature = next(iter(iho_layer))
north_sea_geom = north_sea_feature.GetGeometryRef().Clone()
north_sea_geom.AssignSpatialReference(EPSG4326)
north_sea_envelope = north_sea_geom.GetEnvelope()  # (minX, maxX, minY, maxY)
iho_ds = None  # close datasource; geometry is cloned so it remains valid

print("North Sea envelope (lon_min, lon_max, lat_min, lat_max):", north_sea_envelope)

EPSG4326 = osr.SpatialReference() EPSG4326.ImportFromEPSG(4326) # Force lon/lat axis ordering for consistency with GDAL <3.x and the GeoKit example. EPSG4326.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER) iho_ds = ogr.Open(path_to_sea_basins_shape_file) iho_layer = iho_ds.GetLayer() iho_layer.SetAttributeFilter("NAME = 'North Sea'") north_sea_feature = next(iter(iho_layer)) north_sea_geom = north_sea_feature.GetGeometryRef().Clone() north_sea_geom.AssignSpatialReference(EPSG4326) north_sea_envelope = north_sea_geom.GetEnvelope() # (minX, maxX, minY, maxY) iho_ds = None # close datasource; geometry is cloned so it remains valid print("North Sea envelope (lon_min, lon_max, lat_min, lat_max):", north_sea_envelope)

North Sea envelope (lon_min, lon_max, lat_min, lat_max): (-4.445370674117186, 12.005942344563664, 50.995364665547, 61.01702284842584)

2.2 Display the Region of the North Sea¶

In [4]:

loc-plotting

def _ogr_polygon_to_mpl_patches(geom):
    """Convert an OGR (Multi)Polygon to a list of matplotlib polygon patches."""
    patches = []
    name = geom.GetGeometryName()
    if name == "MULTIPOLYGON":
        polygons = [geom.GetGeometryRef(i) for i in range(geom.GetGeometryCount())]
    else:
        polygons = [geom]
    for poly in polygons:
        # ring 0 = outer, the rest are holes; we ignore holes for simple plots
        ring = poly.GetGeometryRef(0)
        coords = np.array(ring.GetPoints())[:, :2]
        patches.append(MplPolygon(coords, closed=True))
    return patches


fig1, ax1 = plt.subplots(figsize=(6, 6))
ax1.add_collection(
    PatchCollection(_ogr_polygon_to_mpl_patches(north_sea_geom), facecolor="#bcd5f0", edgecolor="black", linewidth=0.5)
)
ax1.set_xlim(north_sea_envelope[0], north_sea_envelope[1])
ax1.set_ylim(north_sea_envelope[2], north_sea_envelope[3])
ax1.set_aspect("equal")
ax1.set_xlabel("Longitude °")
ax1.set_ylabel("Latitude °")
ax1.set_title("North Sea basin (GDAL/OGR)")
plt.show()

def _ogr_polygon_to_mpl_patches(geom): """Convert an OGR (Multi)Polygon to a list of matplotlib polygon patches.""" patches = [] name = geom.GetGeometryName() if name == "MULTIPOLYGON": polygons = [geom.GetGeometryRef(i) for i in range(geom.GetGeometryCount())] else: polygons = [geom] for poly in polygons: # ring 0 = outer, the rest are holes; we ignore holes for simple plots ring = poly.GetGeometryRef(0) coords = np.array(ring.GetPoints())[:, :2] patches.append(MplPolygon(coords, closed=True)) return patches fig1, ax1 = plt.subplots(figsize=(6, 6)) ax1.add_collection( PatchCollection(_ogr_polygon_to_mpl_patches(north_sea_geom), facecolor="#bcd5f0", edgecolor="black", linewidth=0.5) ) ax1.set_xlim(north_sea_envelope[0], north_sea_envelope[1]) ax1.set_ylim(north_sea_envelope[2], north_sea_envelope[3]) ax1.set_aspect("equal") ax1.set_xlabel("Longitude °") ax1.set_ylabel("Latitude °") ax1.set_title("North Sea basin (GDAL/OGR)") plt.show()

3. Match Bathymetry Raster Data to the North Sea¶

3.1 Filter raster tiles overlapping the North Sea¶

GeoKit’s Extent.filterSources returns those tiles whose envelope intersects the extent. We do the same with OGR envelope-intersection logic.

In [5]:

loc-workflow

def envelopes_intersect(a, b):
    # a, b given as (minX, maxX, minY, maxY) — OGR convention
    return not (a[1] < b[0] or a[0] > b[1] or a[3] < b[2] or a[2] > b[3])


all_tile_paths = sorted(pathlib.Path(bathymetry_rasters_dir_str).glob("*.tif"))
bathymetry_datasets = []
for tile_path in all_tile_paths:
    ds = gdal.Open(str(tile_path))
    gt = ds.GetGeoTransform()
    nx, ny = ds.RasterXSize, ds.RasterYSize
    minx = gt[0]
    maxx = gt[0] + nx * gt[1]
    maxy = gt[3]
    miny = gt[3] + ny * gt[5]
    tile_env = (minx, maxx, miny, maxy)
    if envelopes_intersect(tile_env, north_sea_envelope):
        bathymetry_datasets.append(str(tile_path))
    ds = None

bathymetry_datasets

def envelopes_intersect(a, b): # a, b given as (minX, maxX, minY, maxY) — OGR convention return not (a[1] < b[0] or a[0] > b[1] or a[3] < b[2] or a[2] > b[3]) all_tile_paths = sorted(pathlib.Path(bathymetry_rasters_dir_str).glob("*.tif")) bathymetry_datasets = [] for tile_path in all_tile_paths: ds = gdal.Open(str(tile_path)) gt = ds.GetGeoTransform() nx, ny = ds.RasterXSize, ds.RasterYSize minx = gt[0] maxx = gt[0] + nx * gt[1] maxy = gt[3] miny = gt[3] + ny * gt[5] tile_env = (minx, maxx, miny, maxy) if envelopes_intersect(tile_env, north_sea_envelope): bathymetry_datasets.append(str(tile_path)) ds = None bathymetry_datasets

Out[5]:

['geokit/data/gebco_tiles/gebco_2020_n90.0_s-90.0_w-180.0_e180.0_03_08.tif',
 'geokit/data/gebco_tiles/gebco_2020_n90.0_s-90.0_w-180.0_e180.0_03_09.tif',
 'geokit/data/gebco_tiles/gebco_2020_n90.0_s-90.0_w-180.0_e180.0_04_08.tif',
 'geokit/data/gebco_tiles/gebco_2020_n90.0_s-90.0_w-180.0_e180.0_04_09.tif']

3.2 Combine all relevant tiles into a single raster¶

GeoKit’s combineSimilarRasters first verifies that the tiles share SRS / pixel size / aligned origins, then stitches them. The straight-GDAL equivalent is to build a virtual raster (gdal.BuildVRT) and then materialise it as a GeoTIFF (gdal.Translate).

In [6]:

loc-workflow

path_to_bathymetry_file_str = str(current_working_directory.joinpath("bathymetry_raster_NorthSea_gdal.tif"))

vrt_path = str(current_working_directory.joinpath("bathymetry_NorthSea.vrt"))
vrt_ds = gdal.BuildVRT(vrt_path, bathymetry_datasets)
vrt_ds.FlushCache()
vrt_ds = None

merged_ds = gdal.Translate(
    path_to_bathymetry_file_str,
    vrt_path,
    format="GTiff",
    creationOptions=["COMPRESS=LZW", "TILED=YES"],
)
merged_gt = merged_ds.GetGeoTransform()
merged_nx, merged_ny = merged_ds.RasterXSize, merged_ds.RasterYSize
pixel_width = merged_gt[1]
pixel_height = abs(merged_gt[5])
merged_ds = None

assert np.isclose(pixel_width, pixel_height), "Pixel width and height must be equal."
print(f"Merged raster: {merged_nx} x {merged_ny} px, pixel size = {pixel_width} °")

path_to_bathymetry_file_str = str(current_working_directory.joinpath("bathymetry_raster_NorthSea_gdal.tif")) vrt_path = str(current_working_directory.joinpath("bathymetry_NorthSea.vrt")) vrt_ds = gdal.BuildVRT(vrt_path, bathymetry_datasets) vrt_ds.FlushCache() vrt_ds = None merged_ds = gdal.Translate( path_to_bathymetry_file_str, vrt_path, format="GTiff", creationOptions=["COMPRESS=LZW", "TILED=YES"], ) merged_gt = merged_ds.GetGeoTransform() merged_nx, merged_ny = merged_ds.RasterXSize, merged_ds.RasterYSize pixel_width = merged_gt[1] pixel_height = abs(merged_gt[5]) merged_ds = None assert np.isclose(pixel_width, pixel_height), "Pixel width and height must be equal." print(f"Merged raster: {merged_nx} x {merged_ny} px, pixel size = {pixel_width} °")

Merged raster: 10800 x 5400 px, pixel size = 0.004166666666666667 °

3.3 Clip the merged raster to the North Sea polygon¶

GeoKit’s RegionMask.warp aligns and masks the input raster to the polygon. With raw GDAL we use gdal.Warp with cutlineDSName / cropToCutline, telling GDAL to read the cutline geometry from an in-memory OGR datasource.

In [7]:

loc-workflow

# Persist the single North Sea polygon as a GeoJSON cutline file.
# (gdal.Warp cannot resolve in-memory OGR datasources by name — it needs a path.)
cutline_path = str(current_working_directory.joinpath("north_sea_cutline.geojson"))
if pathlib.Path(cutline_path).exists():
    pathlib.Path(cutline_path).unlink()

geojson_drv = ogr.GetDriverByName("GeoJSON")
cutline_ds = geojson_drv.CreateDataSource(cutline_path)
cutline_layer = cutline_ds.CreateLayer("cutline", srs=EPSG4326, geom_type=ogr.wkbPolygon)
feat = ogr.Feature(cutline_layer.GetLayerDefn())
feat.SetGeometry(north_sea_geom.Clone())
cutline_layer.CreateFeature(feat)
feat = None
cutline_ds = None  # flush to disk

warped_path = str(current_working_directory.joinpath("bathymetry_NorthSea_warped_gdal.tif"))
warp_options = gdal.WarpOptions(
    format="GTiff",
    cutlineDSName=cutline_path,
    cropToCutline=True,
    dstSRS="EPSG:4326",
    xRes=pixel_width,
    yRes=pixel_height,
    targetAlignedPixels=True,
    dstNodata=-9999,
    creationOptions=["COMPRESS=LZW", "TILED=YES"],
)
warped_ds = gdal.Warp(warped_path, path_to_bathymetry_file_str, options=warp_options)
warped_ds.FlushCache()

warped_band = warped_ds.GetRasterBand(1)
warped_array = warped_band.ReadAsArray().astype("float64")
warped_nodata = warped_band.GetNoDataValue()
warped_gt = warped_ds.GetGeoTransform()
warped_nx, warped_ny = warped_ds.RasterXSize, warped_ds.RasterYSize
warped_ds = None

if warped_nodata is not None:
    warped_array = np.where(warped_array == warped_nodata, np.nan, warped_array)

print("Warped raster shape:", warped_array.shape)

# Persist the single North Sea polygon as a GeoJSON cutline file. # (gdal.Warp cannot resolve in-memory OGR datasources by name — it needs a path.) cutline_path = str(current_working_directory.joinpath("north_sea_cutline.geojson")) if pathlib.Path(cutline_path).exists(): pathlib.Path(cutline_path).unlink() geojson_drv = ogr.GetDriverByName("GeoJSON") cutline_ds = geojson_drv.CreateDataSource(cutline_path) cutline_layer = cutline_ds.CreateLayer("cutline", srs=EPSG4326, geom_type=ogr.wkbPolygon) feat = ogr.Feature(cutline_layer.GetLayerDefn()) feat.SetGeometry(north_sea_geom.Clone()) cutline_layer.CreateFeature(feat) feat = None cutline_ds = None # flush to disk warped_path = str(current_working_directory.joinpath("bathymetry_NorthSea_warped_gdal.tif")) warp_options = gdal.WarpOptions( format="GTiff", cutlineDSName=cutline_path, cropToCutline=True, dstSRS="EPSG:4326", xRes=pixel_width, yRes=pixel_height, targetAlignedPixels=True, dstNodata=-9999, creationOptions=["COMPRESS=LZW", "TILED=YES"], ) warped_ds = gdal.Warp(warped_path, path_to_bathymetry_file_str, options=warp_options) warped_ds.FlushCache() warped_band = warped_ds.GetRasterBand(1) warped_array = warped_band.ReadAsArray().astype("float64") warped_nodata = warped_band.GetNoDataValue() warped_gt = warped_ds.GetGeoTransform() warped_nx, warped_ny = warped_ds.RasterXSize, warped_ds.RasterYSize warped_ds = None if warped_nodata is not None: warped_array = np.where(warped_array == warped_nodata, np.nan, warped_array) print("Warped raster shape:", warped_array.shape)

Warped raster shape: (2407, 3949)

3.4 Display the bathymetry inside the North Sea¶

We extract the four-corner extent (extent=[xmin, xmax, ymin, ymax]) from the geotransform and use imshow directly. GeoKit’s drawRaster would do the same internally.

In [8]:

loc-plotting

warped_xmin = warped_gt[0]
warped_xmax = warped_gt[0] + warped_nx * warped_gt[1]
warped_ymax = warped_gt[3]
warped_ymin = warped_gt[3] + warped_ny * warped_gt[5]
warped_extent = [warped_xmin, warped_xmax, warped_ymin, warped_ymax]

fig2, ax2 = plt.subplots(figsize=(6, 6))
im = ax2.imshow(
    warped_array,
    cmap="YlGnBu",
    vmin=-600,
    vmax=200,
    extent=warped_extent,
    origin="upper",
)
ax2.add_collection(
    PatchCollection(_ogr_polygon_to_mpl_patches(north_sea_geom), facecolor="none", edgecolor="black", linewidth=0.5)
)
ax2.set_xlim(warped_xmin, warped_xmax)
ax2.set_ylim(warped_ymin, warped_ymax)
ax2.set_aspect("equal")
ax2.set_xlabel("Longitude °")
ax2.set_ylabel("Latitude °")
fig2.colorbar(im, ax=ax2, label="Water Depth (m)")
plt.show()

warped_xmin = warped_gt[0] warped_xmax = warped_gt[0] + warped_nx * warped_gt[1] warped_ymax = warped_gt[3] warped_ymin = warped_gt[3] + warped_ny * warped_gt[5] warped_extent = [warped_xmin, warped_xmax, warped_ymin, warped_ymax] fig2, ax2 = plt.subplots(figsize=(6, 6)) im = ax2.imshow( warped_array, cmap="YlGnBu", vmin=-600, vmax=200, extent=warped_extent, origin="upper", ) ax2.add_collection( PatchCollection(_ogr_polygon_to_mpl_patches(north_sea_geom), facecolor="none", edgecolor="black", linewidth=0.5) ) ax2.set_xlim(warped_xmin, warped_xmax) ax2.set_ylim(warped_ymin, warped_ymax) ax2.set_aspect("equal") ax2.set_xlabel("Longitude °") ax2.set_ylabel("Latitude °") fig2.colorbar(im, ax=ax2, label="Water Depth (m)") plt.show()

4. Load and Display the Offshore Wind Turbine Locations¶

We open the points shapefile with OGR and pull out each feature’s coordinates and FID. GeoKit’s extractFeatures returns a GeoPandas-like dataframe; here we use plain numpy arrays.

4.1 Load Offshore Wind Turbine Locations¶

In [9]:

loc-workflow

turbines_ds = ogr.Open(offshore_wind_points_path_str)
turbines_layer = turbines_ds.GetLayer()

turbine_xy = []
for feat in turbines_layer:
    g = feat.GetGeometryRef()
    turbine_xy.append((g.GetX(), g.GetY()))
turbine_xy = np.array(turbine_xy)
turbines_ds = None

print(f"Loaded {len(turbine_xy)} turbine locations.")

turbines_ds = ogr.Open(offshore_wind_points_path_str) turbines_layer = turbines_ds.GetLayer() turbine_xy = [] for feat in turbines_layer: g = feat.GetGeometryRef() turbine_xy.append((g.GetX(), g.GetY())) turbine_xy = np.array(turbine_xy) turbines_ds = None print(f"Loaded {len(turbine_xy)} turbine locations.")

Loaded 19605 turbine locations.

4.2 Plot Offshore Wind Turbine Locations¶

In [10]:

loc-plotting

fig3, ax3 = plt.subplots(figsize=(6, 6))
im = ax3.imshow(warped_array, cmap="YlGnBu", vmin=-600, vmax=200, extent=warped_extent, origin="upper")
ax3.add_collection(
    PatchCollection(_ogr_polygon_to_mpl_patches(north_sea_geom), facecolor="none", edgecolor="black", linewidth=0.5)
)
ax3.scatter(turbine_xy[:, 0], turbine_xy[:, 1], s=2, c="red", marker=",")
ax3.set_xlim(warped_xmin, warped_xmax)
ax3.set_ylim(warped_ymin, warped_ymax)
ax3.set_aspect("equal")
ax3.set_xlabel("Longitude °")
ax3.set_ylabel("Latitude °")
fig3.colorbar(im, ax=ax3, label="Water Depth (m)")
plt.show()

fig3, ax3 = plt.subplots(figsize=(6, 6)) im = ax3.imshow(warped_array, cmap="YlGnBu", vmin=-600, vmax=200, extent=warped_extent, origin="upper") ax3.add_collection( PatchCollection(_ogr_polygon_to_mpl_patches(north_sea_geom), facecolor="none", edgecolor="black", linewidth=0.5) ) ax3.scatter(turbine_xy[:, 0], turbine_xy[:, 1], s=2, c="red", marker=",") ax3.set_xlim(warped_xmin, warped_xmax) ax3.set_ylim(warped_ymin, warped_ymax) ax3.set_aspect("equal") ax3.set_xlabel("Longitude °") ax3.set_ylabel("Latitude °") fig3.colorbar(im, ax=ax3, label="Water Depth (m)") plt.show()

5. Extract Water Depths at Turbine Locations¶

GeoKit’s interpolateValues(..., mode="near") (the default) returns the nearest pixel value at each point. We do the same with raw GDAL by inverting the geotransform with gdal.InvGeoTransform and indexing into the array we already loaded into memory.

In [11]:

loc-workflow

inv_gt = gdal.InvGeoTransform(warped_gt)

depths = np.full(len(turbine_xy), np.nan)
for i, (lon, lat) in enumerate(turbine_xy):
    px_f, py_f = gdal.ApplyGeoTransform(inv_gt, float(lon), float(lat))
    px, py = int(px_f), int(py_f)
    if 0 <= px < warped_nx and 0 <= py < warped_ny:
        depths[i] = warped_array[py, px]

print("depth statistics (m):")
print(f"  min = {np.nanmin(depths):.1f}")
print(f"  max = {np.nanmax(depths):.1f}")
print(f"  mean= {np.nanmean(depths):.1f}")
depths[:10]

inv_gt = gdal.InvGeoTransform(warped_gt) depths = np.full(len(turbine_xy), np.nan) for i, (lon, lat) in enumerate(turbine_xy): px_f, py_f = gdal.ApplyGeoTransform(inv_gt, float(lon), float(lat)) px, py = int(px_f), int(py_f) if 0 <= px < warped_nx and 0 <= py < warped_ny: depths[i] = warped_array[py, px] print("depth statistics (m):") print(f" min = {np.nanmin(depths):.1f}") print(f" max = {np.nanmax(depths):.1f}") print(f" mean= {np.nanmean(depths):.1f}") depths[:10]

depth statistics (m):
  min = -283.0
  max = -1.0
  mean= -43.0

Out[11]:

array([-32., -32., -30., -27., -25., -24., -26., -28., -26., -24.])

The depths array now holds the same information that the GeoKit notebook stores in the depth column of offshore_wind_points. Combine it with the turbine_xy array to build a table for downstream analysis (cost models, foundation type selection, etc.).