## brcon2026 - 2026-07-31 title: Open analysis and visualization of geospatial data author: Anders Damsgaard (adc) contact: anders@andersdamsgaard.dk gopher://andersdamsgaard.dk https://andersdamsgaard.dk slides: gopher://andersdamsgaard.dk/tmp/brcon2026_adc.md ## About me Work: Geological Survey of Denmark and Greenland (GEUS) Groundwater mapping, geophysics, geological modelling. Part of the work is teaching geophysics to geological surveys and water authorities in African countries. ## Motivation SSC: Strategic Sector Cooperation https://um.dk/en/danida/partners/strategic-sector-cooperation-new #pause SSC is about knowledge and methods exchange. Water: Geophysical instruments produce geospatial data. The instrument costs around EUR 50 000. The software to inspect what it recorded costs a yearly license per seat. #pause The data is public, funded by development aid, and maps drinking water resources. ## Proprietary visualization tools Link a GIS core with tailored 2D and 3D views - Map view (2D) - Depth transects along a line (2D) - 3D views with primitives All require 3D geospatial datarepresentation and handling GIS component is the hard part, compared to the geophysical analysis and visualization tools. Mainly Windows only, bundle database backends, and invent their own GIS. ## GIS troubles: Coordinate reference systems (CRS) The EPSG registry is maintained by the oil industry, originally the European Petroleum Survey Group, now IOGP. EPSG:4326 WGS 84, lat/lon EPSG:32637 WGS 84 / UTM zone 37N (Kenya) EPSG:25832 ETRS89 / UTM zone 32N (Denmark) 11,734 CRSs in PROJ 9.5.1. ## More ways to get coordinates wrong * Axis order. EPSG:4326 is defined as (lat, lon). Most implementations use (lon, lat). GeoJSON mandates lon, lat. CRS84 exists to state longitude first explicitly. #pause * Plate tectonics. ETRS89 is fixed to the Eurasian plate, WGS84/ITRF is not. They diverge by about 2.5 cm/yr, so a 1994 survey and a 2026 survey differ by roughly a metre. ## Elevation The z in POINT Z requires a stated reference surface: ellipsoidal height above a reference ellipsoid orthometric height above the geoid DVR90, EGM2008, ... national and global realizations #pause In Denmark the two differ by roughly 35 to 40 m. Mixing them produces a correctly rendered section at the wrong elevation. #pause There is no single elevation value, only elevation relative to a stated surface. #pause * Post-glacial rebound. Scandinavia is still rising by up to 1 cm/yr. Reference points move vertically. #pause Do not reinvent what others have made, but use FOSS tools. ## Why open formats and tools matter for this work * Groundwater is a drought resilient and under-used source for drinking water in many African countries * TEM is one of the cheaper methods for mapping and understanding water flow in an aquifer system * Widely successful in DK, where drinking water is 99% groundwater * Instruments and training are funded by development programmes * Software licenses are annual, per seat, and priced in EUR #pause Training twenty people requires twenty licenses, or one shared machine. #pause Projects end, licenses expire, and the data remains unused. Knowledge exchange and open formats determine whether the work *will* continue afterwards. ## Transient electromagnetics - TEM method Lay out a wire loop on the ground surface, drive a current through it, then switch the current off abruptly. Observe the induced magnetic signal from the ground. transmitter loop (3 x 3 m or 40 x 40 m) ========================================= ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ surface . o o . t = 10 us . o o . t = 100 us . o o . t = 1 ms . o o . #pause The collapsing field induces current filaments in the ground. They diffuse downward and outward, propagating as smoke rings. Signal intensity reduces over time, and later time represents deeper depths. Resitivity [ohm*m]: a measure of a material's ability to resist the flow of electric current ## The measurement A decay curve. 20-50 float values per sounding, one per time gate. log dB/dt | -6| x | x -8| x x | x x conductive ground: -10| x x slow decay, strong signal | x x -12| x x +-------------------------------------------- -5 -4 -3 -2 log t [s] first gate 12 us last gate 11 ms #pause That is the entire measurement for a place on the earth surface. The rest is inference. The instrument is typically moved along the ground surface or flown through the air to get many measurements. ## Inversion to geology observed decay ---> ??? ---> resistivity vs. depth #pause Assume a layered model, compute the decay it would produce, compare with the observation, adjust, repeat. Gauss-Newton with regularization. #pause The problem is ill-posed: many different models fit the same data equally well. ## Non-uniqueness model A model B ------------------ ------------------ 30 ohm*m, 20 m 30 ohm*m, 20 m ------------------ ------------------ 2 ohm*m, 5 m 1 ohm*m, 10 m both 10 S ------------------ ------------------ 100 ohm*m 100 ohm*m #pause TEM resolves the conductance (thickness / resistivity) of a conductive layer, not thickness and resistivity separately. #pause The regularization determines which of the two models is shown. A plotted model is the data combined with an assumption about smoothness. ## Depth of investigation Below a certain depth the data no longer constrains the model. The inversion still returns numbers. The depth of investigation (DOI) marks where that happens. #pause elev 0 m ### ### ### ### ### ### clay 5 ohm*m -20 ### ### ### ### ### ::: sand 100 ohm*m -40 ::: ::: ::: ::: ::: ... bedrock 4000 ohm*m -60 ::: ::: ::: ... ... -80 ... ... ... ... ... -100 - - DOI - - - - - - - - -120 ? ? ? ? ? s1 s2 s3 s4 s5 #pause The DOI is exported as its own geometry, next to the resitivity result. ## Current state Visualization requires proprietary software by companies, trying to increase profit by selling software that goes along with the instrumentation. This can act as a negative feedback on method application: - against FOSS mindset - closed, centralized development - niche market means high pricing and limited support - out of budget for occasional users or in developing countries #pause The lock-in follows from the business model. It is rarely a deliberate strategy. ## Solution for TEM data GIS is the hard part. Visualization of geophysical data is much simpler. - Transform the data from binary/sql dumps into tabularized clear text - Represent applicable data into geospatial representation - Well-known text (WKT) #pause Store numerical data in flat ASCII files. The Open Geospatial Consortium (OGC) arrived at the same conclusion in 1999. ## Text exports Most vendor software has a text export, usually intended for exchange or debugging. /TEMImage_Beta Models by TEMcompany, version: 4.0.4.6 /Unit Resistivity: [ohm*m], Thickness: [m] /Dummy Resistivity:99999, Dummy Thickness: 99999 /Project Line StationNo InversionName X Y Z DOI \ DataResidual NumLayers Res_001 Res_002 ... Res_030 \ Thick_001 Thick_002 ... Thick_029 #pause /COORDINATE SYSTEM /WGS 84 UTM zone 32N (epsg:32632) /NUMBER OF LAYERS /35 /GATE TIMES (s) / 1.2450E-5 1.5510E-5 1.9320E-5 ... 1.1130E-2 / LINE_NO UTMX UTMY DATE TIME RECORD ELEVATION \ NUMDATA SEGMENTS RESDATA RESTOTAL RHO_1 ... RHO_35 ## The XYZ format Two vendors, three export dialects, the same physics. * Comment lines start with '/', in most of the variants * Missing data is 9999, 99999, or -9999 * The data model is encoded in column names: Res_001..Res_030 * 35 layers gives 105 columns * Header keys on one line, values on the next * Fixed-width, except where it is whitespace delimited * CRS sometimes given, sometimes not. #pause The wide layout comes from spreadsheet exports. Parsing it is straightforward once the variants are known. ## Well-known text format: Geospatial primitives POINT (x y) POINT Z (x y z) LINESTRING (x y, x y, ...) POLYGON ((x y, ..., x y)) first == last MULTIPOINT ((x y), (x y)) MULTILINESTRING ((x y, x y), (...)) MULTIPOLYGON (((x y, ...)), ((...))) GEOMETRYCOLLECTION (POINT (x y), ...) #pause That is the complete geometry model. OGC 06-103r4 / ISO 19125. https://docs.ogc.org/is/06-103r4/06-103r4.pdf Anything shown on a map is one of these eight types. ## How the data becomes geospatial primitives Geometry is written as WKT strings into the output and parsed by QGIS on load. For a sounding at (x, y) with surface elevation z and depth of investigation doi: void print2dpoint(FILE *fp, double x, double y) { fprintf(fp, "POINT (%f %f)", x, y); } #pause void print3dpoint(FILE *fp, double x, double y, double z) { fprintf(fp, "POINT Z (%f %f %f)", x, y, z); } ## How the data becomes geospatial primitives (cont.) Each resistivity layer is a vertical line segment between the top and bottom elevation of that layer: void print3dline(FILE *fp, double x1, double y1, double z1, double x2, double y2, double z2) { fprintf(fp, "LINESTRING Z (%f %f %f, %f %f %f)", x1, y1, z1, x2, y2, z2); } #pause z_top = z - depth_top z_bot = z - depth_bottom Coordinates within a point are space separated. Points are comma separated. Mixing the two separators is a common error. ## Fence diagrams One sounding is a stack of LINESTRING Z segments, coloured by resistivity. Placed along a survey line, the stacks form a section. #pause z ^ ------------ Ground ------------- | ### ### ### ### ### ### ### | ### ### ### ### ### ### ### | ::: ::: ### ### ### ::: ::: | ::: ::: ::: ::: ::: ::: ::: | ... ::: ::: ::: ::: ::: ... | ... ... ... ... ... ... ... +----------------------------------------> distance ## Without the plugin $ { printf 'wkt,res\n' awk 'NR>4 {printf "\"POINT Z (%s %s %s)\",%s\n", \ $5, $6, $7, $11}' stem.xyz } > soundings.csv The same WKT string is understood by: PostGIS ST_GeomFromText('POINT Z (...)') SpatiaLite same GDAL/OGR CSV driver, GeoJSON driver, and others shapely, JTS, GEOS, and most web map libraries ## GeoJSON WKT is the text grammar. GeoJSON is the same model in JSON. POINT Z (543210 9950000 1180) #pause { "type": "Feature", "geometry": { "type": "Point", "coordinates": [543210, 9950000, 1180] }, "properties": {"line": "L01", "rho": 12.4} } #pause RFC 7946 (2016): The GeoJSON format. Round-trips losslessly with WKT. The RFC mandates a single coordinate reference system: WGS 84, longitude before latitude. This removes a common source of errors. ## WKT, GeoJSON, or GeoPackage WKT compact, greppable, no attributes of its own GeoJSON self-describing, nests, widely supported GeoPackage sqlite(3), indexed, queryable, binary, one file ## Rasters Everything so far has been vector data: geometry with coordinates. The other half of geospatial data is raster: a grid of values with a known position. vector raster ------ ------ POINT Z (x y z) [ 12 14 13 11 ] LINESTRING Z (...) [ 14 18 17 12 ] POLYGON ((...)) [ 15 19 22 16 ] where is the feature [ 13 15 14 12 ] what is the value here #pause Terrain models, satellite imagery, geophysical grids, model output, scanned maps. #pause The plugin can correct sounding elevations by sampling band 1 of a DEM raster, so the section sits on the measured terrain. ## GeoTIFF TIFF dates from 1986. It is a container of tagged records, and the tag space is extensible. In 1994 Niles Ritter at NASA JPL placed the coordinate system and georeferencing into a set of private TIFF tags. 33550 ModelPixelScale 33922 ModelTiepoint 34264 ModelTransformation 34735 GeoKeyDirectory the CRS is defined here #pause The result is a georeferenced image that is still a valid TIFF. $ xdg-open dem.tif # opens as an ordinary image $ gdalinfo dem.tif # also reports the location #pause GeoTIFF spread because existing TIFF readers kept working. ## The geotransform An affine transform from pixel coordinates to world coordinates: x = gt[0] + col*gt[1] + row*gt[2] y = gt[3] + col*gt[4] + row*gt[5] gt[0], gt[3] origin (upper left corner) gt[1], gt[5] pixel size (gt[5] is usually negative) gt[2], gt[4] rotation, usually 0 #pause $ gdalinfo dem.tif | grep -A1 'Origin\|Pixel Size' #pause Six doubles and a CRS definition are enough to georeference a raster. The remainder of the format is compression and tiling. ## Cloud Optimized GeoTIFF A COG is a GeoTIFF with three additional requirements: 1. internally tiled rather than stripped 2. overviews (downsampled pyramids) included 3. a documented layout, with the header first #pause $ gdal_translate in.tif out.tif -of COG \ -co COMPRESS=ZSTD -co PREDICTOR=2 #pause This then works against a static file on an ordinary HTTP server: $ gdalinfo /vsicurl/https://host/huge.tif $ gdal_translate -projwin ... /vsicurl/https://host/huge.tif \ small.tif #pause A 200 GB file, a 2 MB transfer, using HTTP range requests. No server software beyond a static file server, and no API. ## VRT A VRT file describes a raster mosaic in XML. EPSG:25832 ... tile_0001.tif ... $ gdalbuildvrt denmark.vrt tiles/*.tif $ ls -la denmark.vrt -rw-r--r-- 1 adc adc 742431 denmark.vrt #pause Many tiles and TBs on disk, described by one XML file that can be read and kept in version control. - COG: makes one raster efficient to read, especially over HTTP. - VRT: makes many rasters appear as one dataset. ## Raster and vector conversion GDAL ships the conversions as separate programs: gdal_rasterize vector -> raster gdal_polygonize raster -> vector gdal_contour raster -> LINESTRING contours gdaldem hillshade DEM -> shaded relief gdallocationinfo raster -> value at a coordinate #pause # terrain elevation at a sounding position $ gdallocationinfo -valonly -geoloc dem.tif 543210 9950000 ## The open geospatial stack PROJ 1983, coordinate transformations, C GEOS port of Java Topology Suite (JTS), geometry predicates, C++ GDAL/OGR format translation, C++ libtiff GeoTIFF is a TIFF SQLite GeoPackage is a SQLite database QGIS 2002, Gary Sherman, C++/Qt, GPL #pause Commercial GIS products generally link against the same libraries. The differentiation is mostly in the user interface. ## ogr2ogr # inspect a dataset $ ogrinfo -so -al data.gpkg # convert between formats $ ogr2ogr -f GeoJSON out.json in.shp # reproject during conversion $ ogr2ogr -t_srs EPSG:4326 wgs84.gpkg utm.gpkg # filter with SQL, without a database server $ ogr2ogr -sql "SELECT * FROM layers WHERE rho < 10" \ clay.gpkg soundings.gpkg ## File formats and archival Shapefile (ESRI, 1998, dBASE III attribute table from 1983) - 2 GB limit per file - field names limited to 10 characters - no NULL values, no proper date/time type - one shapefile is 3 to 8 separate files - still the default export in most software #pause GeoPackage (OGC, 2014) - one file, SQLite, documented, no size limit - readable with sqlite3(1) #pause CSV with WKT - readable with cat(1), diff(1) and git(1) GeoTIFF for raster data - a 1986 TIFF still decodes today ## QGIS GIS is the hard part.... -> Use an FOSS GIS platform and expand it for the geophysical data analysis and visualization. QGIS: https://github.com/qgis/QGIS First release in 2002 GPL licensed. QGIS 4 release: * Qt 6 replaces Qt 5 * deprecated APIs ease plugin migration * few new user-facing features The map, the formats and the command line remain the same: GDAL/OGR + PROJ + GEOS + GeoPackage ## The QGIS plugin for TEM data in SSC tem_loader/core.py 590 lines no QGIS imports tem_loader/interpolate.py 380 lines IDW interpolation tem_loader/tem_loader.py 879 lines Qt dialogs test/ 2600 lines #pause * core.py has no QGIS dependency, so it can be unit tested with pytest and called from a shell script * no numpy, scipy, or pandas * three vendor dialects normalized into one internal schema * output is a GeoPackage written next to the source file #pause https://plugins.qgis.org/plugins/tem_loader/ Around 1000 installs so far around the world. Used every day in Denmark, Kenya, Marrocco, Ethiopia, Ukraine, France, ... ## UNIX principles applied Small is beautiful 590 lines, one parser Make each program do one thing parse, style, render Build a prototype immediately printf("POINT (%f %f)") Choose portability text Store numerical data in ASCII WKT Use software leverage GDAL, PROJ and GEOS solve the hard parts already Avoid captive user interfaces the plugin is optional Make every program a filter xyz2wkt < in > out ## Summary * Vendor lock-in in science is largely a file format problem * Writing the parser is easy * WKT and GeoJSON cover the vector side, GeoTIFF the raster side * GDAL, PROJ, GEOS, QGIS and SQLite are already available * Coordinate reference systems remain the main source of error Blog post with some screenshots: https://adamsgaard.dk/tem-loader.html Plugin page on QGIS.org: https://plugins.qgis.org/plugins/tem_loader/ Source code: git://src.adamsgaard.dk/qgis-tem-loader ## Thanks for listening