Genereren van 3D Tiles: i3dm

Laatst bijgewerkt: 10 maart 2025

• download van BGT de objecten als GML bestand.
  https://app.pdok.nl/lv/bgt/download-viewer/
  
• laad het GML bestand in een Postgres database met ogr2ogr.
  ga naar het ogr2ogr programma
  
  ogr2ogr -update -append -f “PostgreSQL” PG:”host=localhost port=5432 dbname=i3dm user=postgres password=***** active_schema=public” “D:\bgt_nederland\bgt_vegetatieobject.gml” -progress -t_srs EPSG:28992 –config CPL_DEBUG ON
  
  Opmerkingen
  – omdat we werken met spatial data moet je in Postgres de Postgis extensie installeren
  – ogr2ogr is een opensource tool waarmee data geconverteerd kan worden tussen verschillende bronnen.
  ogr2ogr maakt onderdeel uit van GDAL en kan zelfstandig gedownload en geïnstalleerd worden. De GDAL library wordt ook automatisch meegeleverd met QGIS.
  
• de srid is verloren gegaan tijdens de import
  herstel deze als volgt:
  
  UPDATE solitaryvegetationobject
  SET geometrie2dvegetatieobject = ST_SetSRID(geometrie2dvegetatieobject, 28992)
  ;
  
• We moeten de z-waarde berekenen. Dat doen we in QGIS met behulp van
  een python script.
  Genereer met onderstaande sql statement de input voor dit python script.
  
  create table qgis_input_file
  as
  select ‘QgsPointXY(‘|| ST_X(geometrie2dvegetatieobject)||’, ‘||ST_Y(geometrie2dvegetatieobject)||’),’
  from solitaryvegetationobject
  WHERE ST_GeometryType(geometrie2dvegetatieobject) = ‘ST_Point’
  –limit 100
  ;
  
  Je kunt de uitkomst kopiëren en plakken in het python script. Wanneer er sprake is van grote hoeveelheden
  dan kun je met ogr2ogr een csv dump maken.
  
  ogr2ogr -f “CSV” \LAP48007\Users\xxx\QGIS_input_file.csv PG:”host=localhost port=5432 dbname=i3dm user=postgres password=***** active_schema=public” -sql “SELECT * FROM qgis_input_file” -progress
  
• Het python script
  
  #1. Rasterlaag ophalen
  layer = iface.activeLayer()
  provider = layer.dataProvider()
  
  2. Lijst met X, Y coördinaten (in RD-coördinaten)#
  punten = [
  QgsPointXY(48112.164017905976, 417624.9924212946),
  QgsPointXY(44036.66601821994, 417729.7163859736),
  QgsPointXY(44067.63901837748, 417784.47238490544),
  QgsPointXY(44120.230017926646, 417808.3303892127)
  ]
  
  #3. Open een bestand om resultaten op te slaan
  output_file = ‘C:/Users/path/output.sql’
  with open(output_file, “w”) as file:
  
  # 4. Loop door alle punten en haal hoogtewaarden op
  for point in punten:
  z_value = provider.sample(point, 1)[0] # Alleen de hoogte (Z)
  
  line = f”INSERT INTO qgis_import (geom, scale, rotation, model) ” \
  f”VALUES (ST_GeomFromText(‘POINT({point.x()} {point.y()} {z_value})’, 28992), 10, 10, ‘tree.glb’);\n”
  
  #print(line.strip())
  file.write(line)
  print(f”Hoogtewaarden opgeslagen in: {output_file}”)
  
• Het python script voert een raster scan uit en bepaalt voor elk coördinaat de Z-waarde.
  Dit gebeurt met behulp van de AHN geotiff bestanden. (Algemeen Hoogtebestand Nederland)
  
  Deze bestanden kunnen gedownload worden vanaf PDOK
  
  – Open vervolgens QGIS
  – Navigeer: Kaartlagen > Laag toevoegen > Rasterlaag
  – Open een python console: CTRL+Alt+P
  – Klik op het icoon ‘Editor Tonen’
  – Open of plak het python script
  – Voer het script uit
  
  
• Maak een nieuwe tabel aan in Postgres
  
  CREATE TABLE IF NOT EXISTS public.qgis_import
  (
  id integer NOT NULL DEFAULT nextval(‘mydata_id_seq’::regclass),
  geom geometry(PointZ,28992),
  scale double precision,
  scale_non_uniform double precision[],
  rotation double precision,
  model character varying COLLATE pg_catalog.”default”,
  tags json,
  CONSTRAINT qgis_import_pkey PRIMARY KEY (id)
  ;
  
• Het resultaat zal er als volgt uitzien. Het import sql statement bevat X, Y en Z.
  
  INSERT INTO qgis_import (geom, scale, rotation, model) VALUES (ST_GeomFromText(‘POINT(108756.658 513904.316 0.08529999852180481)’, 28992), 10, 10, ‘tree.glb’);
  INSERT INTO qgis_import (geom, scale, rotation, model) VALUES (ST_GeomFromText(‘POINT(108755.762 513904.519 0.12210000306367874)’, 28992), 10, 10, ‘tree.glb’);
  INSERT INTO qgis_import (geom, scale, rotation, model) VALUES (ST_GeomFromText(‘POINT(108820.742 513850.342 0.06159999966621399)’, 28992), 10, 10, ‘tree.glb’);
  
  Dit bovenstaande bestand kan weer ingeladen worden in Postgres.
  
• De bovenstaande insert statement zijn gebaseerd op het projectiestelsel EPSG:28992.
  Cesium verwacht echter EPSG:4326.
  Daarom moeten we nog een transformatie uitvoeren.
  We doen dit wederom in een nieuwe tabel.
  
  create table tiling_table
  as
  select id
  , ST_transform (geom, 4326) geom
  , scale
  , scale_non_uniform
  , rotation, model
  , tags
  from qgis_import
  ;
  
  Voorzie deze tabel van een index. Dit verbetert de performance enorm.
  
  CREATE INDEX geom_idx
  ON tiling_table
  USING GIST (geom);
  
• De tabel ’tiling_table’ vormt de input voor het programma i3dm.export
  
  dotnet run — -c “Host=localhost;Username=postgres;Password=*****;Database=BAG;Port=5432” -t tiling_table -o “C:\Users\afdmu\www\3dtilesnederland\tiles”