When dealing with geo data, there are essentially two major types of data: vector data and raster data. Vector data includes points, lines and polygons, while raster data is represented as a two- or three-dimensional pixel map.
If your aim is to map every fast-food restaurant in the US, you would probably represent each of them as a point, i.e. you will deal with vector data.
Raster data, on the other hand, could, for example, be the geographical altitude of every single pixel within some area. Or, in the case of Orbio – where we largely work with satellite imagery – a satellite observation of a certain geographical area is also represented as raster of pixels, each of which holds, for instance, a value for the observed light intensity / brightness.
Every raster has a certain spatial resolution, that is, every pixel represents a certain area in m². Moreover, raster can – but do not have to – be geo-referenced, i.e. aligned to some geographical coordinate system (the most prominent of which is EPSG:4326).
PostGIS, Geoserver and WMS
The by far most popular way to store and query any kind of geo data is using PostGIS, which is a brilliant, open-source extension to the PostgreSQL database management system. It supports vector data types (e.g. POINT, POLYGON, LINESTRING) as well as raster data (RASTER type).
This is where WMS (aka. web map service) comes in. WMS is a protocol, that specifies how to request geo data on the web. Input parameters include a bounding box of the requested area, an SRID code and others, while the output of a WMS service is rendered images (usually PNGs, JPEGs or TIFFs).
Consequently, as the developer of a GIS application, you somehow need to provide a WMS endpoint, that, on one end, fetches data from PostGIS, and, on the other end, delivers it to the web frontend.
One approach to do so is to use a map server, like Geoserver, Mapserver or QGIS Server. Geoserver is arguably the most feature-rich software among those, however, unfortunately, it inherently lacks the ability to serve raster data from a PostGIS data base.
Custom WMS endpoint implementation
Since WMS itself is not that complex as a protocol, you can quite easily come up with your own WMS-compatible endpoint (at least one that suffices the GetMap capability).
PL/pgSQL database routine
The first component is a quite comprehensive SQL function, that does most of the heavy lifting. Once created in your database it can be called as part of SQL queries with a couple of parameters. Connect to your PostGIS database and run the following script to initially create the function.
EXECUTE 'SELECT ST_Buffer($1, 20);' USING var_env INTO var_env_buf;
var_sql := 'SELECT srid, ST_AsRaster($4,$5,$6,pixel_types,nodata_values,nodata_values) As erast FROM raster_columns WHERE r_table_schema = $1AND r_table_name = $2AND r_raster_column=$3'; EXECUTE var_sql INTO var_srid, var_erast USING param_schema, param_table, 'rast', var_env, param_width, param_height; var_sql := 'WITH r AS (SELECT ST_Clip(rast,' || CASE WHEN var_srid = param_srid THEN '$7' ELSE 'ST_Transform($7,$2)' END || ') As rast FROM' || quote_ident(param_schema) || '.' || quote_ident(param_table) || ' WHERE ST_Intersects(rast,' || CASE WHEN var_srid = param_srid THEN '$3' ELSE 'ST_Transform($3,$2)' END || ') limit15) SELECT ST_Clip(ST_Union(rast), $3) As rast FROM (SELECT ST_Resample(' || CASE WHEN var_srid = param_srid THEN 'rast' ELSE 'ST_Transform(rast,$1)' END || ',$6,true,''NearestNeighbor'') As rast FROM r) Asfinal'; EXECUTE var_sql INTO var_result USING param_srid, var_srid, var_env, param_width, param_height, var_erast, var_env_buf; var_sql := 'SELECT ST_MapAlgebra($1, $2, ''[rast2]'', ''8BUI''::text, ''FIRST'', ''[rast2]'', NULL::text) rast'; EXECUTE var_sql INTO var_result USING var_erast, var_result; IF var_result IS NULL THEN var_result := var_erast; END IF; RETURN CASE WHEN param_format ILIKE 'image/jpeg' THEN ST_AsJPEG(ST_ColorMap(var_result, 1, 'bluered')) ELSE ST_AsPNG(ST_ColorMap(var_result, 1, 'bluered')) END; END; $$;
The function takes a set of parameters, which mostly correspond to WMS parameters and it returns a fully rendered PNG (or JPEG) image as a byte array.
Second step is to implement the actual web endpoint, that can be called from GIS clients using HTTP. The following examples show implementations of a WMS endpoint (utilizing the above get_rast_tile SQL function) in two of the most popular Python web frameworks, FastAPI and Flask.
Both examples assume the presence of a PostGIS database called geodb, which includes a table observations for storing raster data. The table has to, at least, consist of a primary key column and a column rast of type RASTER.
You can either create that table manually (not recommended), or use the command-line program raster2pgsql automatically create an appropriate table structure.
import uvicorn as uvicorn from fastapi import FastAPI, Query, Depends from sqlalchemy import text from sqlalchemy.ext.asyncio import create_async_engine, AsyncSession from sqlalchemy.orm import sessionmaker from starlette.responses import StreamingResponse
# Service method for retrieving WMS data from database defget_wms(db: Session, params: typing.Dict[str, typing.Any]) -> typing.BinaryIO: with db.begin(): params['crs'] = int(params['crs'].removeprefix('EPSG:'))
# Required for ST_AsPNG and ST_AsJPEG to work db.execute("SET postgis.gdal_enabled_drivers = 'ENABLE_ALL';") db.execute("SET postgis.enable_outdb_rasters TO true;")