Geospatial Data Storage, Interpolation & GIS Export for Environmental IoT #

By the time environmental telemetry reaches this stage it has already been ingested, timestamp-aligned, calibrated, and quality-flagged — but it is still just a stream of points: a device ID, a coordinate, a timestamp, and a value. The people who actually make decisions with that data — air quality regulators, hydrologists, urban planners, ecologists — do not work in Kafka topics or pandas frames. They work in QGIS, they overlay layers on basemaps, they draw contours and compute zonal statistics, and they expect a continuous surface, not a scatter of dots. The engineering problem this section solves is the translation: taking clean point telemetry and turning it into stored, queryable, interpolated, and shareable geospatial products that a GIS analyst can open without a Python interpreter anywhere in sight.

That translation has four moving parts — a spatial store that indexes millions of readings for fast bounding-box and nearest-neighbor queries, an interpolation layer that estimates a continuous field from sparse sensors, a vector export path for points and contours, and a raster export path for gridded surfaces. Each has its own failure modes, and the CRS discipline that binds them together is where most pipelines quietly go wrong.

Pipeline Architecture: From Point Telemetry to Geospatial Products #

The diagram below shows the horizontal flow this section builds. Calibrated point telemetry lands in a spatially indexed PostGIS store; queries pull a point set for a target time; an interpolation stage estimates a continuous surface; and two export paths deliver a vector product for GIS desktops and a cloud-optimized raster for tiling and analysis.

The one invariant that runs through all four stages is the coordinate reference system. Storage keeps a canonical WGS84 copy for portability and a projected copy for computation; interpolation must run in a metric projected CRS; vector export must land back in EPSG:4326; raster export must carry an explicit transform and CRS tag. Lose track of the CRS at any stage and the whole product is silently wrong. The sections below take each stage in turn.

Core Concept A — Spatial Storage and Indexing #

The store is the foundation. Environmental networks accumulate readings relentlessly — a few hundred nodes at one-minute cadence is tens of millions of rows a year — and every downstream query is spatial: “all PM2.5 readings inside this city boundary in the last hour,” “the ten nearest active sensors to this point,” “the bounding box covering catchment X.” Without a spatial index those queries degrade into full table scans that make interactive GIS impossible.

PostGIS is the default answer. It stores geometry natively, supports GiST and BRIN indexes over that geometry, and speaks the same SQL your ingestion layer already writes. The PostGIS storage and spatial indexing for sensor networks approach covers schema design, index selection, and bulk loading in depth; the minimal pattern is a table with a typed geometry column and a spatial index over it:

# python 3.11 · psycopg2-binary==2.9.9 · geopandas==0.14.3 · shapely==2.0.4
# pyproj==3.6.1 · scipy==1.12.0 · pykrige==1.7.1 · xarray==2024.3.0
# rioxarray==0.15.0 · rasterio==1.3.9
import psycopg2

DDL = """
CREATE TABLE IF NOT EXISTS sensor_reading (
    id          BIGINT GENERATED ALWAYS AS IDENTITY PRIMARY KEY,
    device_id   TEXT        NOT NULL,
    reading_ts  TIMESTAMPTZ NOT NULL,
    variable    TEXT        NOT NULL,          -- 'pm25', 'temperature', ...
    value       DOUBLE PRECISION,
    qc_flag     SMALLINT    NOT NULL DEFAULT 0,-- 0=ok, propagated from QC stage
    geom        GEOMETRY(Point, 4326) NOT NULL -- canonical WGS84 lon/lat
);
CREATE INDEX IF NOT EXISTS sensor_reading_geom_gist
    ON sensor_reading USING GIST (geom);
CREATE INDEX IF NOT EXISTS sensor_reading_ts_brin
    ON sensor_reading USING BRIN (reading_ts);
"""

with psycopg2.connect("dbname=env_iot") as conn, conn.cursor() as cur:
    cur.execute("CREATE EXTENSION IF NOT EXISTS postgis;")
    cur.execute(DDL)

The index choice is deliberate. A GiST index over geom accelerates the && bounding-box operator and ST_DWithin proximity queries that dominate spatial lookups. A BRIN index over reading_ts is a near-free win for append-only time-ordered data: because rows arrive in roughly ascending timestamp order, BRIN summarizes each block with a min/max range and prunes entire blocks for time-window queries at a fraction of a B-tree’s storage cost. The trade-off between GiST and BRIN for the geometry column itself — and when append order justifies BRIN over geometry — is exactly the decision covered under GiST versus BRIN indexes for sensor geometry.

The qc_flag column is not decoration. It carries forward the verdict from the automated calibration, validation and anomaly detection stage so that interpolation can filter to trusted readings without a join back to a separate QC table. Store the flag next to the value; interpolate only rows where it is clean.

Core Concept B — Spatial Interpolation #

Point telemetry answers “what was the value at this sensor.” A GIS analyst almost always wants “what was the value everywhere” — a continuous field they can contour, threshold, or overlay. Interpolation estimates that field from the sparse points, and the two workhorses are inverse distance weighting and kriging. The spatial interpolation with kriging and IDW for sensor fields treatment works through variogram fitting, neighbor selection, and validation; the essential mechanics of IDW over a target grid are short:

# numpy==1.26.4 · scipy==1.12.0 · pyproj==3.6.1
import numpy as np
from scipy.spatial import cKDTree
from pyproj import Transformer

def idw_surface(lonlat: np.ndarray, values: np.ndarray,
                grid_lon: np.ndarray, grid_lat: np.ndarray,
                utm_epsg: int, power: float = 2.0, k: int = 12) -> np.ndarray:
    """IDW interpolation performed in a projected (metric) CRS.

    lonlat: (n, 2) WGS84 lon/lat of sensors. values: (n,) readings.
    grid_lon/grid_lat: 2-D meshgrid of output cell centres in WGS84.
    Returns a 2-D array of interpolated values on the grid.
    """
    fwd = Transformer.from_crs(4326, utm_epsg, always_xy=True)
    sx, sy = fwd.transform(lonlat[:, 0], lonlat[:, 1])          # sensors -> metres
    gx, gy = fwd.transform(grid_lon.ravel(), grid_lat.ravel())  # grid -> metres

    tree = cKDTree(np.column_stack([sx, sy]))
    dist, idx = tree.query(np.column_stack([gx, gy]), k=k)
    dist = np.where(dist == 0, 1e-9, dist)                      # avoid divide-by-zero
    w = 1.0 / dist ** power
    out = np.sum(w * values[idx], axis=1) / np.sum(w, axis=1)
    return out.reshape(grid_lon.shape)

Note the transform to UTM before any distance is computed. Distances in raw longitude/latitude degrees are meaningless as weights — a degree of longitude is ~111 km at the equator and ~55 km at 60° latitude — so interpolating on degrees stretches the surface east-west by latitude. Project to metres, interpolate, and the weights reflect true ground distance.

Kriging goes further: it fits a variogram model to the empirical spatial correlation and produces both an estimate and a per-cell variance, which is invaluable for showing analysts where the surface is trustworthy and where it is guesswork. That extra rigor costs fitting time and demands enough well-distributed points, which is the crux of choosing ordinary kriging versus IDW for sparse sensor networks. For a low-latency air quality heatmap refreshed every few minutes, building an IDW air quality surface from point sensors is usually the pragmatic starting point. Interpolation is also the natural bridge from the real-time stream processing and spatial analytics layer, where a live surface may be recomputed on each window close.

Core Concept C — Vector and Raster Export #

An interpolated surface trapped in a NumPy array helps no one. The last stage serializes it into the two formats GIS tooling actually consumes: GeoJSON for vector features and GeoTIFF for rasters. These are not interchangeable — points, sensor metadata, and contour lines belong in vector; the continuous grid belongs in raster — and a complete product usually ships both.

The vector path is where GeoJSON and QGIS export workflows live. GeoJSON is a strict format: WGS84 only, longitude before latitude. The GeoJSON and QGIS export workflows cover styling, property schemas, and contour generation; the core is reprojecting to EPSG:4326 and writing with the GeoJSON driver:

# geopandas==0.14.3 · shapely==2.0.4
import geopandas as gpd
from shapely.geometry import Point

def export_points_geojson(df, utm_epsg: int, path: str) -> None:
    """Serialize QC-clean sensor readings to WGS84 GeoJSON for QGIS."""
    gdf = gpd.GeoDataFrame(
        df.loc[df["qc_flag"] == 0].copy(),
        geometry=[Point(xy) for xy in zip(df["lon"], df["lat"])],
        crs="EPSG:4326",
    )
    # Guarantee output CRS regardless of how the frame was built.
    gdf = gdf.to_crs(4326)
    gdf.to_file(path, driver="GeoJSON")

Carrying the qc_flag and provenance columns into the GeoJSON properties lets analysts filter and symbolize in QGIS directly — the pattern detailed in exporting QC-flagged data to GeoJSON for QGIS.

The raster path serializes the interpolated grid. Wrapping the array in an xarray.DataArray with named x/y coordinates and handing it to rioxarray gives you a georeferenced GeoTIFF with a single call, and the gridded raster fields with xarray and rioxarray approach explains dimension ordering, nodata handling, and cloud-optimized output:

# xarray==2024.3.0 · rioxarray==0.15.0 · rasterio==1.3.9 · numpy==1.26.4
import numpy as np
import xarray as xr
import rioxarray  # noqa: F401 — registers the .rio accessor

def grid_to_geotiff(surface: np.ndarray, lon: np.ndarray, lat: np.ndarray,
                    path: str, nodata: float = -9999.0) -> None:
    """Write a 2-D interpolated field to a georeferenced GeoTIFF.

    lat MUST be descending (north-to-south) so the affine transform
    matches GDAL's top-left origin convention.
    """
    da = xr.DataArray(
        np.where(np.isnan(surface), nodata, surface),
        dims=("y", "x"),
        coords={"y": lat, "x": lon},
    )
    da.rio.write_crs("EPSG:4326", inplace=True)
    da.rio.write_nodata(nodata, inplace=True)
    da.rio.to_raster(path, driver="GTiff", tiled=True, compress="deflate")

Setting tiled=True plus internal overviews is what turns an ordinary GeoTIFF into a cloud-optimized one that QGIS and web tile servers can pan without reading the whole file — the subject of writing cloud-optimized GeoTIFFs from sensor grids.

Production Implementation Patterns #

One canonical CRS, plus a working projection #

Adopt a two-CRS discipline and enforce it everywhere. Store the canonical geometry in EPSG:4326 — it is portable, it is what GeoJSON demands, and it is unambiguous across tools. Do every metric computation (interpolation, buffering, distance, area) in a projected CRS, typically the local UTM zone. The failure this prevents is subtle: code that “works” on degrees produces plausible-looking but distorted surfaces that no gate catches. Reproject explicitly at each boundary and never let a bare (lon, lat) pair enter a distance calculation.

Filter on QC flag before interpolation, not after #

The clean-data contract from the upstream QC stage exists to be used. Interpolating raw readings and masking the surface afterward is backwards — a single anomalous spike will pull every neighboring interpolated cell toward it before you ever mask anything. Filter qc_flag == 0 at the point-set query, so the interpolator only ever sees trusted inputs. This is also why the flag is stored inline in PostGIS rather than in a side table.

Preserve provenance through every serialization #

Each exported feature and each raster should carry enough metadata to reconstruct how it was made: the variable name, the units, the timestamp or window, the interpolation method and its parameters, the source CRS, and the sensor count that fed the surface. In GeoJSON this goes in feature properties or a top-level metadata object; in GeoTIFF it goes in tags via da.rio.update_tags(). Analysts who open a raster six months later have no other way to know whether they are looking at an IDW power-2 surface from eight sensors or a kriged field from forty.

Vectorize the grid, never loop cells #

Interpolation output grids are large — a 500×500 field is 250,000 cells. Build the target grid with numpy.meshgrid, transform coordinate arrays in a single pyproj call, and let the cKDTree query resolve all grid points at once. Per-cell Python loops turn a sub-second operation into minutes and are the most common performance regression in interpolation code.

Operationalizing and Scaling #

Recompute on a schedule, or on demand #

Surfaces are derived products, and you have two ways to keep them current. Precompute on a fixed cadence — regenerate the citywide PM2.5 raster every five minutes from the latest window — when a fixed set of standard products is served to many viewers; the cost is bounded and cache-friendly. Compute on demand — interpolate only when an analyst requests a specific variable, time, and bounding box — when the parameter space is large and most combinations are never requested. Most deployments do both: precompute the handful of dashboard surfaces, and expose an on-demand endpoint for ad-hoc analysis that reads the point set straight from the indexed PostGIS store.

Partition and retain the point store #

A single sensor_reading table grows without bound. Partition it by time (monthly range partitions) so that queries for a recent window touch only recent partitions and old data can be detached and archived to object storage as Parquet in one operation. The BRIN index on reading_ts complements partitioning rather than replacing it. Keep the raw calibrated points as the system of record and treat every raster and GeoJSON as a reproducible derivative that can be regenerated and therefore safely expired.

Version the derived products #

When a calibration coefficient is revised upstream, or the interpolation parameters change, previously exported surfaces become inconsistent with newly generated ones. Tag each exported product with the code version, the interpolation parameters, and the input point-set hash. Storing derivatives in object storage under a content-addressed or run-tagged prefix lets you reproduce exactly which points and settings produced a given raster — the same reproducibility discipline the real-time stream processing and spatial analytics pipelines apply to their windowed outputs. This is what makes an air quality surface defensible in a regulatory review months later.

Failure Modes and Gotchas #

Interpolating in geographic degrees #

The single most damaging mistake in this pipeline is running IDW or kriging on raw longitude/latitude. Distances in degrees are not isotropic — a degree east covers far less ground than a degree north at high latitude — so a surface interpolated on degrees is compressed north-south relative to reality, and every distance-weighted or variogram-based estimate is wrong in a way that looks completely plausible. Always project to a metric CRS first.

Extrapolation beyond the sensor hull #

IDW and kriging estimate between sensors. Ask them for values far outside the convex hull of the network — the corners of a rectangular output grid that extends past the outermost sensors — and they return confident numbers backed by nothing. Clip the output grid to the sensor hull buffered by a sensible radius, or mask cells whose nearest sensor exceeds a maximum distance, so the product never implies coverage the network does not have.

NaN and nodata confusion in rasters #

A grid cell with no valid estimate must be written as an explicit nodata value, and that value must be declared on the raster. Leave NaNs in the array without calling write_nodata, and different readers disagree: some render NaN as black, some as white, some propagate it through zonal statistics and poison the result. Choose a sentinel outside the physical range of the variable (never 0 for a quantity that can legitimately be 0) and tag it.

Axis order and flipped rasters #

GeoJSON is longitude-then-latitude; several projected formats and some libraries are the reverse. Combined with a y-coordinate written ascending instead of descending, this produces rasters that load upside down or vector layers that plot in the wrong hemisphere near (0, 0). Assert always_xy=True on every pyproj transformer, confirm the raster’s latitude coordinate descends, and spot-check one known sensor location in QGIS before trusting a batch.

Spatial autocorrelation assumptions near point sources #

Interpolation assumes nearby locations have similar values. Near a strong point source — a smokestack, a busy intersection — the gradient is steep and that assumption fails; the interpolated value between a high sensor and a low one represents neither. Where the network straddles known point sources, prefer denser sampling and nearest-neighbor estimation over smooth interpolation, and treat any smooth surface across such a gradient as indicative rather than accurate.

Frequently Asked Questions #

Should I store sensor coordinates as geography or geometry in PostGIS?

Store as geometry in a projected CRS when your network fits inside one UTM zone — distance and buffer operations are faster and interpolation happens in projected space anyway. Use geography only when readings span continents and you need true great-circle distances without picking a projection. For most environmental deployments, geometry in the local UTM zone plus a WGS84 copy for export is the pragmatic choice.

When should I use kriging instead of inverse distance weighting?

Use ordinary kriging when you need per-cell uncertainty estimates and your field shows a fittable spatial correlation structure — typical of slowly varying quantities like temperature or groundwater level. Use IDW when you need low-latency surfaces, have dense coverage, or the field is dominated by local point sources where variogram fitting is unstable. Kriging is slower and requires enough well-distributed points to fit a reliable variogram.

Which CRS should the interpolation run in?

Run interpolation in a projected CRS whose units are metres — a local UTM zone or a national grid — so that distance-based weights and variogram lags are computed in real ground distance. Interpolating directly on WGS84 degrees distorts distances by latitude and produces stretched surfaces. Reproject points into the projected CRS, interpolate, then reproject the output grid back to EPSG:4326 for GeoJSON export.

How do I make a raster from sensor readings open correctly in QGIS?

Write the grid as a GeoTIFF with an explicit CRS, a correct affine transform, and a defined nodata value using rioxarray. The most common reason a raster loads with no georeference or upside down is a missing or flipped transform — ensure the y-coordinate is written in descending order or set the transform explicitly. Add internal tiling and overviews to produce a cloud-optimized GeoTIFF that QGIS can pan smoothly.

Why does my exported GeoJSON have coordinates in the wrong place?

GeoJSON is defined in WGS84 longitude/latitude (EPSG:4326) with longitude first. If features appear off the coast of Africa near (0, 0) or land in the wrong hemisphere, you exported projected metre coordinates without reprojecting, or you swapped lat and lon. Always reproject the GeoDataFrame to EPSG:4326 before calling to_file with the GeoJSON driver, and confirm the axis order is lon, lat.

---

Related #

• PostGIS Storage & Spatial Indexing for Sensor Networks — schema design, GiST and BRIN indexes, and bulk-loading millions of readings
• Spatial Interpolation with Kriging & IDW for Sensor Fields — variogram fitting, neighbor selection, and choosing a method for sparse networks
• GeoJSON & QGIS Export Workflows — WGS84 serialization, property schemas, contours, and QGIS styling
• Gridded Raster Fields with xarray & rioxarray — DataArray grids, nodata handling, and cloud-optimized GeoTIFF output
• IoT Sensor Data Ingestion & Spatial Synchronization — upstream stage: MQTT and Kafka ingestion, CRS mapping, and timestamp alignment
• Automated Calibration, Validation & Anomaly Detection — upstream QC that supplies the calibrated, quality-flagged points this stage stores
• Real-Time Stream Processing & Spatial Analytics — sibling stage: windowed aggregation and live surface recomputation on stream events