{"id":2304,"date":"2023-09-01T16:56:06","date_gmt":"2023-09-01T14:56:06","guid":{"rendered":"https:\/\/www.speich.net\/articles\/?p=2304"},"modified":"2023-11-13T10:44:47","modified_gmt":"2023-11-13T09:44:47","slug":"vector-tiles-with-postgis-php-and-openlayers","status":"publish","type":"post","link":"https:\/\/www.speich.net\/articles\/en\/2023\/09\/01\/vector-tiles-with-postgis-php-and-openlayers\/","title":{"rendered":"Vector tiles with PostGIS, PHP and OpenLayers"},"content":{"rendered":"\n

This is an example about how to create vector tiles with PostGIS 3.3<\/a> (and PostgreSQL 15.3), serve them with PHP 8.1<\/a> and use them in OpenLayers 7.5<\/a>. In contrast to the few tutorials available on the internet, it not only uses the latest PostGIS version 3.3, which further simplified the generation of vector tiles by introducing the new function ST_TileEnvelope<\/a>, but also shows the use case, when you want to create the tiles in a projection other than the default Spherical Mercator coordinate system (EPSG:3857<\/a>).<\/p>\n\n\n\n

\"Screenshot
The production regions of Switzerland rendered with OpenLayers. The vector tiles are created on the fly with PostGIS and php. The clipped tile boundaries are shown on purpose.<\/figcaption><\/figure>\n\n\n\n\n\n\n\n

Create the tiles with a SQL query<\/h2>\n\n\n\n
WITH\n    \/* Prepare your data, in my case:\n            - ST_Union: aggregate multiple polygons into one per id\n            - ST_Transform: transform geometries to the desired target projection srid:21781\n            - ST_ForcePolygonCCW: make sure they are rotated counter clock wise\n     *\/\n    mydata AS (\n        SELECT ST_ForcePolygonCCW(\n            ST_Transform(ST_Union(t.geom), 21781)\n        ) geom, someColumn\n        FROM someTable t\n        GROUP by t.someColumn\n    ),\n\n    \/* Create the tile envelope with ST_TileEnvelope for use in the function ST_AsMVTGeom below.\n     * If you want to use your own projection (EPSG:21781), create the bounds with ST_MakeEnvelope.\n     * Remember, to use bounds which are preferably square instead of rectangular.\n     *\/\n    bounds AS (\n        SELECT ST_TileEnvelope(:z, :x, :y, ST_MakeEnvelope(485072, -52106, 837120, 299942, 21781)) AS geom\n    ),\n\n    \/*\n     * Transform the geometry with ST_AsMVTGeom into the coordinate space of a Mapbox Vector Tile (MVT)\n     * clipping it to tile bounds if necessary.\n     * Note: for the sake of this example we set the buffer to zero, to show the clipped boundaries\n     *\/\n    mvtgeom AS (\n        SELECT ST_AsMVTGeom(mydata.geom, bounds.geom, 4096, 0) AS geom,\n            otherColumn\n        FROM mydata, bounds\n        WHERE ST_Intersects(mydata.geom, bounds.geom)\n    )\n\/*\n * Create the binary MVT format with ST_AsMVT. All row columns other than the geometry are encoded as feature attributes.\n *\/\nSELECT ST_AsMVT(t.*, 'myLayerName') mvt\nFROM mvtgeom t<\/code><\/pre>\n\n\n\n
If you don’t need to use a projection other than Web Mercator, then the above SQL query simplifies to:<\/p>\n\n\n\n
WITH mydata AS (\n        SELECT ST_ForcePolygonCCW(\n            ST_Transform(ST_Union(t.geom), 3857)\n        ) geom, prodreg otherColumn\n        FROM someTable t\n        GROUP by t.otherColumn\n    ),\n    bounds AS (\n        SELECT ST_TileEnvelope(:z, :x, :y) AS geom\n    ),\n    mvtgeom AS (\n        SELECT ST_AsMVTGeom(mydata.geom, bounds.geom) AS geom,\n               otherColumn\n        FROM mydata, bounds\n        WHERE ST_Intersects(data.geom, bounds.geom)\n    )\nSELECT ST_AsMVT(t.*, 'myLayerName') mvt\nFROM mvtgeom t<\/code><\/pre>\n\n\n\n
Serve the tiles with PHP<\/h2>\n\n\n\n
Now, to access the PostgreSQL database, we use the database abstraction layer PHP Data Objects<\/a> (PDO). The database driver specific to PostgreSQL has be installed separately, on (Ubuntu) Linux with:<\/p>\n\n\n\n
$ apt install php8.1-pgsql<\/code><\/p>\n\n\n\n
For this example, we create the resource \/my-webroot\/service.php<\/p>\n\n\n\n
\/\/ service.php accessed as \/myWebroot\/service.php\/{z}\/{x}\/{y}.pbf\n\n\/\/ Grab the tile indices from the path:\n[$z, $x, $y] = explode('\/', $_SERVER['PATH_INFO'])\n$y = substr($y, 0, -4);\n\n\/\/ Query the tiles with the above SQL:\n$sql = \"WITH data ... FROM mvtgeom t\";\n$db = connect();\n$stmt = $db->prepare($sql);\n\n\/\/ Bind the placeholders to the query:\n$stmt->execute(['x' => $x, 'y' => $y, 'z' => $z]);\n$row = $stmt->fetch();\n\n\/\/ Set the correct content type of the response header:\nheader('Content-Type: application\/vnd.mapbox-vector-tile; utf-8');\n\n\/\/ Output the requested tile\necho stream_get_contents($row['mvt']);<\/code><\/pre>\n\n\n\n
It’s that simple !<\/p>\n\n\n\n
Have the browser cache your tiles<\/h2>\n\n\n\n
To speed things up, query the database to create the tiles only, when you actually have to send a response body, i.e. wrap the database query in a function or method and only call it when the Etag has changed:<\/p>\n\n\n\n
\/\/ Create an Etag header to let the browser cache the tiles (optional):\n\/\/ Note: The identifier needs to be wrapped in double quotes.\n$etag = md5($tileDateLastUpdate)\n$etagHeader = isset($_SERVER['HTTP_IF_NONE_MATCH']) ? trim($_SERVER['HTTP_IF_NONE_MATCH']) : false;\nheader('Etag: \"'.$etag.'\"');\n\/\/ There is no need for the browser to send a request for this tile again until the next day (the response is fresh for 60*60*24 seconds):\nheader('Cache-Control: max-age=86400');\n\n\/\/ If we have a matching Etag, we know that the content was not modified:\nif ($etagHeader === $etag) {    \n    \/\/ Do not query the database here, since we don't send a body at all.\n    http_response_code(304);\n} else {\n    \/\/ Only connect to your database and query it here, since this takes longer:\n    echo myTileService.get($x, $y, $z);\n}<\/code><\/pre>\n\n\n\n
Use the tiles as a vector layer in OpenLayers<\/h2>\n\n\n\n
OpenLayers also uses Spherical Mercator (EPSG:3857) as the default. Since we are using the Swiss CH1903\/LV03 coordinate system in this example (EPSG:21781<\/a>), we have to explicitly set the projection on the view and the tile source:<\/p>\n\n\n\n
import MVT from '.\/node_modules\/ol\/format\/MVT.js';\nimport VectorTileLayer from '.\/node_modules\/ol\/layer\/VectorTile.js';\nimport VectorTileSource from '.\/node_modules\/ol\/source\/VectorTile.js';\nimport View from '.\/node_modules\/ol\/View.js';\nimport Map from '.\/node_modules\/ol\/Map.js';\nimport proj4 from '.\/node_modules\/proj4';\nimport {register} from '.\/node_modules\/ol\/proj\/proj4.js';\n\nconst proj = {\n  \/**\n   * projection definition\n   * CH1903 \/ LV03 -- Swiss CH1903 \/ LV03\n   * @see https:\/\/epsg.io\/21781\n   * @type {String}\n   *\/\n  definition: '+proj=somerc +lat_0=46.9524055555556 +lon_0=7.43958333333333 +k_0=1 +x_0=600000 +y_0=200000 +ellps=bessel +towgs84=674.374,15.056,405.346,0,0,0,0 +units=m +no_defs +type=crs',\n\n  name: 'EPSG:21781',\n  extent: [485071.58, 74261.72, 837119.8, 299941.79],\n  center: [660013.54, 185171.98],\n};\n\nproj4.defs(proj.name, proj.definition);\nregister(proj4);\n\nconst map = new Map({\n  target: 'map-container',\n  view: new View({\n    projection: proj.name,\n    center: proj.center,\n    extent: proj.extent,\n    showFullExtent: true,\n    zoom: 1,\n  }),\n  layers: [\n    new VectorTileLayer({\n      source: new VectorTileSource({\n        projection: proj.name,\n        extent: proj.extent,\n        format: new MVT(),\n        url: 'https:\/\/result-maps.test\/service.php\/tile\/{z}\/{x}\/{y}.pbf',\n        maxZoom: 14,\n        tileSize: 256,\n      })\n    })\n  ]\n});<\/code><\/pre>\n\n\n\n
If we wanted to stick to Spherical Mercator and display the vector tiles over OpenStreetMap, the JavaScript code would simplify to:<\/p>\n\n\n\n
import MVT from '.\/node_modules\/ol\/format\/MVT.js';\nimport VectorTileLayer from '.\/node_modules\/ol\/layer\/VectorTile.js';\nimport VectorTileSource from '.\/node_modules\/ol\/source\/VectorTile.js';\nimport View from '.\/node_modules\/ol\/View.js';\nimport Map from '.\/node_modules\/ol\/Map.js';\nimport TileLayer from '.\/node_modules\/ol\/layer\/Tile.js';\nimport {OSM} from '.\/node_modules\/ol\/source.js';\n\nconst map = new Map({\n  target: 'map-container',\n  view: new View({\n    center: [915602.8072104596, 5911929.4687797595],\n    extent: [663464, 5751550, 1167741, 6075303],\n    showFullExtent: true,\n    zoom: 3,\n  }),\n  layers: [\n    new TileLayer({source: new OSM()}),\n    new VectorTileLayer({\n      source: new VectorTileSource({\n        format: new MVT(),\n        url: 'https:\/\/result-maps.test\/service.php\/tile\/{z}\/{x}\/{y}.pbf',\n        maxZoom: 14,\n      })\n    })\n  ]\n});<\/code><\/pre>\n\n\n\n
Style your vector tiles<\/h2>\n\n\n\n
...\nimport Style from '.\/node_modules\/ol\/style\/Style.js';\nimport Fill from '.\/node_modules\/ol\/style\/Fill.js';\nimport Stroke from '.\/node_modules\/ol\/style\/Stroke.js';\n\nnew VectorTileLayer({\n      source: new VectorTileSource(...),\n\n      \/\/ additionally, style the tiles with a white stroke and an orange, semitransparent fill\n      style: new Style({\n        stroke: new Stroke({\n          color: '#fff',\n          width: 2\n        }),\n        fill: new Fill({\n          color: 'rgba(255, 153, 0, 0.5)'\n        })\n      })\n    })<\/code><\/pre>\n\n\n\n
\"Screenshot
The production regions of Switzerland (polygons) styled with a white stroke and semitransparent orange fill. OpenStreetMap is set as an additional (background) layer.<\/figcaption><\/figure>\n\n\n\n
We can also use the feature properties to color our polygons, in my case we use the property regionid. Instead of passing a style instance to the style property, we use a function:<\/p>\n\n\n\n
    new VectorTileLayer({\n      source: new VectorTileSource(...),\n      \/\/ create a fill color from the feature property region id\n      style: function(feature) {\n        let color;\n        const prop = feature.getProperties();\n        switch (prop.regionid) {\n          case 1:\n            color = 'rgba(213,223,229, 0.6)';\n            break;\n          case 2:\n            color = 'rgba(201,177,189, 0.6)';\n            break;\n          case 3:\n            color = 'rgba(180,149,148, 0.6)';\n            break;\n          case 4:\n            color = 'rgba(127,145,114, 0.6)';\n            break;\n          case 5:\n            color = 'rgba(86,117,104, 0.6)';\n            break;\n        }\n        style.getFill().setColor(color);\n\n        return style;\n      },\n    })<\/code><\/pre>\n\n\n\n
\"Screenshot
The polygons (production regions) can be styled dynamically with a function where the tile features are passed as an argument.<\/figcaption><\/figure>\n\n\n\n
Worth reading<\/h2>\n\n\n\n