# Cartographic application with PostGIS and pg_tileserv
This workshop and related exercices are distributed under a CC-BY-SA licence ![](https://mirrors.creativecommons.org/presskit/buttons/80x15/png/by-sa.png).
Workshop content :
[TOC]
## General presentation
### Goals of this workshop
You will learn the mechanisms to build a cartographic application around PostgreSQL/PostGIS.
You will also learn a few basic French words and frenglish ;-)
### Work environment
We will work on the following environment :
- Pre-installed Virtualbox virtual machine
- Linux Lubuntu 20.04 system
The virtual machine is available here :
- https://share.oslandia.net/public/9a2ceb
- sha1sum : a483d999a6a777324273534e55e7a17271d88554
On the system you will find the following tools pre-installed :
- PostgreSQL 12
- PostGIS 3
- QGIS 3.20
As for the training platform, we use the platform setup by FOSS4G organizers for video-conferencing, screen sharing and chat.
Do not hesitate to :
- Ask any question on the chat
- Ask to go slower
- Ask to explain an item again
- Ask for a break
Should you wish to have the keyboard in your own locale, you can run the following in a terminal (here to french, adapt to your locale name ):
```bash
setxkbmap fr
```
### Credentials
- Ubuntu : oslandia / oslandia
- PostgreSQL : postgres / oslandia (change this with `sudo -u postgres psql -c "alter user postgres password 'oslandia'"`)
Obviously, do not use this virtual machine for production...
### Components
- PostGIS : the spatial database, a PostgreSQL extension
- QGIS : The famous desktop GIS
- pg_tileserv : Serving vector tiles as MVT
- pg_featureserv : Serving vector data as json/geojson
# Workshop steps
Part 1 - data into PostGIS
- Download cadastral parcel data
- Create a PostGIS database
- Integrate the data
Part 2 - pg_tileserv
- Download and install pg_tileserv
- Configure pg_tileserv
- Check if it works
Part 3 - a web application
- Use OpenLayers as a client to access vector tiles
Part 4 - PostGIS for cartographic data processing
- Use PostGIS and pg_tileserv with data processing
## Part 1 - Data to PostGIS
Download the data :
- [Cadastral parcel data for Paris (75)](https://cadastre.data.gouv.fr/data/etalab-cadastre/2020-07-01/shp/departements/75/)
```
mkdir data
cd data
wget -e robots=off -np --mirror https://cadastre.data.gouv.fr/data/etalab-cadastre/2020-07-01/shp/departements/75/
cd cadastre.data.gouv.fr/data/etalab-cadastre/2020-07-01/shp/departements/75
for i in *zip;do unzip $i;done
```
Create a database and activate PostGIS extension.
Create the database :
```shell=sh
sudo -u postgres createdb cadastre
psql -h localhost -U postgres -d cadastre
```
Create the extension :
```sql
CREATE EXTENSION postgis;
```
Load the data into the database :
```shell=sh
# Check what is inside the data file
ogrinfo batiments.shp -so -al
# Load it into the database
SHAPE_ENCODING="LATIN1" ogr2ogr -f "PostgreSQL" PG:"host=localhost user=postgres dbname=cadastre password=oslandia" -nlt PROMOTE_TO_MULTI ./batiments.shp "batiments"
```
Repeat for all other `.shp` files :
- **communes.shp**
- **parcelles.shp**
- **batiments.shp**
- feuilles.shp
- prefixes sections.shp
- sections.shp
- subdivisions_fiscales.shp
### Check loaded data
- Open QGIS
- Add a PostGIS data source
- name : cadastre
- service : leave empty
- host : localhost
- port : 5432
- database : cadastre
- auth : basic
- login : postgres
- password : oslandia ( :warning: you should have setup the database password before, see part "credentials" above)
- Open the buildings table
- Zoom to the data
- Check the attribute table
## Part 2 - Pg_tileserv
### Installation
https://access.crunchydata.com/documentation/pg_tileserv/latest/installation/
```bash
cd
mkdir tileserv
cd tileserv
wget https://postgisftw.s3.amazonaws.com/pg_tileserv_latest_linux.zip
unzip pg_tileserv_latest_linux.zip
```
### run the server to check the installation
```bash
export DATABASE_URL=postgresql://postgres:oslandia@localhost/cadastre
./pg_tileserv
```
Open a browser at the following URL :
http://localhost:7800
In the console, you can see the HTTP queries in the server logs when you visualize the data in browser.
What is the format used for transfer ?
### configuration
We can create a configuration file for pg_tileserv and use it.
Edit `tileserv_cadastre.toml` :
```bash
# Database connection
DbConnection = "user=postgres password=oslandia host=localhost dbname=cadastre"
# Close pooled connections after this interval
DbPoolMaxConnLifeTime = "1h"
# Hold no more than this number of connections in the database pool
DbPoolMaxConns = 4
# Look to read html templates from this directory
AssetsPath = "./assets"
# Accept connections on this subnet (default accepts on all subnets)
HttpHost = "0.0.0.0"
# Accept connections on this port
HttpPort = 7800
# Advertise URLs relative to this server name
# default is to look this up from incoming request headers
# UrlBase = "http://yourserver.com/"
# Resolution to quantize vector tiles to
DefaultResolution = 4096
# Rendering buffer to add to vector tiles
DefaultBuffer = 256
# Limit number of features requested (-1 = no limit)
MaxFeaturesPerTile = -1
# Advertise this minimum zoom level
DefaultMinZoom = 10
# Advertise this maximum zoom level
DefaultMaxZoom = 22
# Allow any page to consume these tiles
CORSOrigins = "*"
# logging information?
Debug = false
```
Restart the server :
```
./pg_tileserv --config tileserv_cadastre.toml
```
## Part 3 - pg_featureserv
### Installation
https://access.crunchydata.com/documentation/pg_featureserv/latest/installation/installing/
```bash
cd
mkdir featureserv
cd featureserv
wget https://postgisftw.s3.amazonaws.com/pg_featureserv_latest_linux.zip
unzip pg_featureserv_latest_linux.zip
```
### run the server to check the installation
```bash
export DATABASE_URL=postgresql://postgres:oslandia@localhost/cadastre
./pg_featureserv
```
Open a browser at the following URL :
http://localhost:9000
In the console, you can see the HTTP queries in the server logs when you visualize the data in browser.
What is the format used for transfer ?
### configuration
We can create a configuration file for pg_featureserv and use it.
:::info
A full example is available in the `config` folder of pg_featureserv archive.
:::
Edit `featureserv_cadastre.toml` :
```bash
[Server]
# Accept connections on this subnet (default accepts on all)
HttpHost = "0.0.0.0"
# IP ports to listen on
HttpPort = 9000
# String to return for Access-Control-Allow-Origin header
# CORSOrigins = "*"
# set Debug to true to run in debug mode (can also be set on cmd-line)
Debug = true
[Database]
# Database connection
# postgresql://username:password@host/dbname
# DATABASE_URL environment variable takes precendence if set.
DbConnection = "postgresql://postgres:oslandia@localhost/cadastre"
[Paging]
# The default number of features in a response
LimitDefault = 20
# Maxium number of features in a response
LimitMax = 10000
[Metadata]
# Title for this service
#Title = "pg-featureserv"
# Description of this service
#Description = "Crunchy Data Feature Server for PostGIS"
```
Restart the server :
```
./pg_featureserv --config featureserv_cadastre.toml
```
And check everything is still working as espected.
## Part 4 - Web application
We will use OpenLayers to visualize our vector tiles.
### Our OpenLayers application
```
cd
mkdir app
cd app
featherpad index.html
```
Edit the `index.html` file and insert the following code :
```htmlmixed=
<!doctype html>
<html lang="en">
<head>
<meta charset="utf-8">
<meta name="viewport" content="width=device-width, initial-scale=1">
<title>OpenLayers Vector Tiles</title>
<!-- CSS/JS for OpenLayers map -->
<link rel="stylesheet" href="https://openlayers.org/en/v6.1.1/css/ol.css" type="text/css">
<script src="https://openlayers.org/en/v6.1.1/build/ol.js"></script>
<style>
body {
padding: 0;
margin: 0;
}
html, body, #map {
height: 100%;
width: 100%;
font-family: sans-serif;
}
#meta {
background-color: rgba(255,255,255,0.75);
color: black;
z-index: 2;
position: absolute;
top: 10px;
left: 20px;
padding: 10px 20px;
margin: 0;
}
</style>
</head>
<body>
<div id="meta">
<h2>Our cartographic application !</h2>
<ul>
<li><a href="http://localhost:7800/"> pg_tileserv Server</a></li>
</ul>
</div>
<div id="map"></div>
<script>
var vectorServer = "http://localhost:7800/";
var vectorSourceLayer = "public.parcelles";
var vectorProps = "?properties=ogc_fid,id,commune,prefixe,section,numero,contenance,created,updated"
var vectorUrl = vectorServer + vectorSourceLayer + "/{z}/{x}/{y}.pbf" + vectorProps;
var vectorStyle = new ol.style.Style({
stroke: new ol.style.Stroke({
width: 2,
color: "#ff00ff99"
}),
fill: new ol.style.Fill({
color: "#ff00ff33"
})
});
var vectorLayer = new ol.layer.VectorTile({
source: new ol.source.VectorTile({
format: new ol.format.MVT(),
url: vectorUrl
}),
style: vectorStyle
});
var baseLayer = new ol.layer.Tile({
source: new ol.source.XYZ({
url: "http://osm.oslandia.io/styles/osm-bright/{z}/{x}/{y}.png"
})
});
var map = new ol.Map({
target: 'map',
view: new ol.View({
center: [260369,6249082],
zoom: 16
}),
layers: [baseLayer, vectorLayer]
});
</script>
</body>
</html>
```
### Run the application
- Start pg_tileserv as above
- Start the OpenLayers application thanks to an HTTP Python micro-server
```shell
cd app
python3 -m http.server
```
If an alternative port is needed for the python http server, it can be specified using the *bind* flag, e.g.
```shell
python3 -m http.server --bind 127.0.0.1 9999
```
- Now open [the application in a browser](http://localhost:8000)
We can now analyze how our application works :
- Analyze network queries with the web developer toolbox in the browser
- changer the rendering style for vector tyles
- add a layer for communes ( Paris is divided into separate communes )
- add a layer for batiments ( = buildings )
```htmlmixed=
<!doctype html>
<html lang="en">
<head>
<meta charset="utf-8">
<meta name="viewport" content="width=device-width, initial-scale=1">
<title>OpenLayers Vector Tiles</title>
<!-- CSS/JS for OpenLayers map -->
<link rel="stylesheet" href="https://openlayers.org/en/v6.1.1/css/ol.css" type="text/css">
<script src="https://openlayers.org/en/v6.1.1/build/ol.js"></script>
<style>
body {
padding: 0;
margin: 0;
}
html, body, #map {
height: 100%;
width: 100%;
font-family: sans-serif;
}
#meta {
background-color: rgba(255,255,255,0.75);
color: black;
z-index: 2;
position: absolute;
top: 10px;
left: 20px;
padding: 10px 20px;
margin: 0;
}
</style>
</head>
<body>
<div id="meta">
<h2>Our cartographic application !</h2>
<ul>
<li><a href="http://localhost:7800/"> pg_tileserv Server</a></li>
<li><a href="http://localhost:9000/"> pg_featurserv Server</a></li>
</ul>
</div>
<div id="map"></div>
<script>
var vectorServer = "http://localhost:7800/";
var parcellesSourceLayer = "public.parcelles";
var parcellesProps = "?properties=ogc_fid,id,commune,prefixe,section,numero,contenance,created,updated"
var communesSourceLayer = "public.communes";
var buildingsSourceLayer = "public.batiments";
var parcellesUrl = vectorServer + parcellesSourceLayer + "/{z}/{x}/{y}.pbf" + parcellesProps;
const communesUrl = vectorServer + communesSourceLayer + "/{z}/{x}/{y}.pbf";
const buildingsUrl = vectorServer + buildingsSourceLayer + "/{z}/{x}/{y}.pbf";
var parcellesStyle = new ol.style.Style({
stroke: new ol.style.Stroke({
width: 1.5,
color: "#907a4f"
}),
fill: new ol.style.Fill({
color: "#e2c07c"
})
});
var buildingsStyle = new ol.style.Style({
stroke: new ol.style.Stroke({
width: 1.5,
color: "#417143"
}),
fill: new ol.style.Fill({
color: "#73c77688"
})
});
// parcelles
var parcellesLayer = new ol.layer.VectorTile({
source: new ol.source.VectorTile({
format: new ol.format.MVT(),
url: parcellesUrl
}),
style: parcellesStyle
});
// communes
var communesLayer = new ol.layer.VectorTile({
source: new ol.source.VectorTile({
format: new ol.format.MVT(),
url: communesUrl
}),
//style: communesStyle
});
// parcelles
var buildingsLayer = new ol.layer.VectorTile({
source: new ol.source.VectorTile({
format: new ol.format.MVT(),
url: buildingsUrl
}),
style: buildingsStyle
});
var baseLayer = new ol.layer.Tile({
source: new ol.source.XYZ({
url: "http://osm.oslandia.io/styles/osm-bright/{z}/{x}/{y}.png"
})
});
var map = new ol.Map({
target: 'map',
view: new ol.View({
center: [260369,6249082],
zoom: 16
}),
layers: [baseLayer, communesLayer, parcellesLayer, buildingsLayer]
});
</script>
</body>
</html>
```
### Our work environment
Now that we have a few components running, it will be interesting to have the following in your work environment
- a terminal with 3 tabs :
- pg_tileserv running ( `./pg_tileserv -c tileserv_cadastre.toml` )
- pg_featureserv running ( `./pg_featureserv -c featureserv_cadastre.toml` )
- python HTTP server running ( `python3 -m http.server` )
- an open SQL console ( `psql -U postgres -d cadastre -h localhost` )
- QGIS Desktop opened with PostGIS database access
- A browser with 2 tabs open :
- pg_tileserv homepage ( http://localhost:7800 )
- our web application ( http://localhost:8000 )
- A file manager open
- A text editor with our web application source files open
### Add a layer switcher
We will now improve our application with a layer switcher widget. Edit the file `index2.html` :
```htmlmixed=
<!doctype html>
<html lang="en">
<head>
<meta charset="utf-8">
<meta name="viewport" content="width=device-width, initial-scale=1">
<title>OpenLayers Vector Tiles</title>
<!-- CSS/JS for OpenLayers map -->
<link rel="stylesheet" href="https://openlayers.org/en/v6.1.1/css/ol.css" type="text/css">
<script src="https://openlayers.org/en/v6.1.1/build/ol.js"></script>
<script src="https://unpkg.com/ol-layerswitcher@3.7.0"></script>
<style>
body {
padding: 0;
margin: 0;
}
html, body, #map {
height: 100%;
width: 100%;
font-family: sans-serif;
}
#meta {
background-color: rgba(255,255,255,0.75);
color: black;
z-index: 2;
position: absolute;
top: 10px;
left: 20px;
padding: 10px 20px;
margin: 0;
}
</style>
<link rel="stylesheet" href="https://unpkg.com/ol-layerswitcher@3.7.0/src/ol-layerswitcher.css" />
</head>
<body>
<div id="meta">
<h2>Our cartographic application !</h2>
<ul>
<li><a href="https://localhost:7800">pg_tileserv server</a></li>
</ul>
</div>
<div id="map"></div>
<script>
var vectorServer = "http://localhost:7800/";
var vectorSourceLayer = "public.parcelles";
var vectorProps = "?properties=ogc_fid,id,commune,prefixe,section,numero,contenance,created,updated"
var vectorUrl = vectorServer + vectorSourceLayer + "/{z}/{x}/{y}.pbf" + vectorProps;
var vectorStyle = new ol.style.Style({
stroke: new ol.style.Stroke({
width: 2,
color: "#ff00ff99"
}),
fill: new ol.style.Fill({
color: "#ff00ff33"
})
});
var vectorLayer = new ol.layer.VectorTile({
title:'Parcelles'
,source: new ol.source.VectorTile({
format: new ol.format.MVT(),
url: vectorUrl
}),
style: vectorStyle
});
var baseLayer = new ol.layer.Tile({
title: 'osm',
source: new ol.source.XYZ({
url: "http://osm.oslandia.io/styles/osm-bright/{z}/{x}/{y}.png"
})
});
var map = new ol.Map({
target: 'map',
view: new ol.View({
center: [260369,6249082],
zoom: 16
}),
layers: [baseLayer, vectorLayer]
});
var layerSwitcher = new ol.control.LayerSwitcher({
tipLabel: 'Legend', // Optional label for button
groupSelectStyle: 'children' // Can be 'children' [default], 'group' or 'none'
});
map.addControl(layerSwitcher);
</script>
</body>
</html>
```
### Add a popup display
We improve our application again, adding a popup when we click on the map. Edit file `index3.html` :
```htmlmixed=
<!doctype html>
<html lang="en">
<head>
<meta charset="utf-8">
<meta name="viewport" content="width=device-width, initial-scale=1">
<title>OpenLayers Vector Tiles</title>
<!-- CSS/JS for OpenLayers map -->
<link rel="stylesheet" href="https://openlayers.org/en/v6.1.1/css/ol.css" type="text/css">
<script src="https://openlayers.org/en/v6.1.1/build/ol.js"></script>
<script src="https://unpkg.com/ol-layerswitcher@3.7.0"></script>
<script src="https://unpkg.com/ol-popup@4.0.0"></script>
<style>
body {
padding: 0;
margin: 0;
}
html, body, #map {
height: 100%;
width: 100%;
font-family: sans-serif;
}
#meta {
background-color: rgba(255,255,255,0.75);
color: black;
z-index: 2;
position: absolute;
top: 10px;
left: 20px;
padding: 10px 20px;
margin: 0;
}
</style>
<link rel="stylesheet" href="https://unpkg.com/ol-layerswitcher@3.7.0/src/ol-layerswitcher.css" />
<link rel="stylesheet" href="https://unpkg.com/ol-popup@4.0.0/src/ol-popup.css" />
</head>
<body>
<div id="meta">
<h2>Our cartographic application !</h2>
<ul>
<li><a href="https://localhost:7800"> pg_tileserv server</a></li>
</ul>
</div>
<div id="map"></div>
<script>
var vectorServer = "http://localhost:7800/";
var vectorSourceLayer = "public.parcelles";
var vectorProps = "?properties=ogc_fid,id,commune,prefixe,section,numero,contenance,created,updated"
var vectorUrl = vectorServer + vectorSourceLayer + "/{z}/{x}/{y}.pbf" + vectorProps;
var vectorStyle = new ol.style.Style({
stroke: new ol.style.Stroke({
width: 2,
color: "#ff00ff99"
}),
fill: new ol.style.Fill({
color: "#ff00ff33"
})
});
var vectorLayer = new ol.layer.VectorTile({
title:'Parcelles'
,source: new ol.source.VectorTile({
format: new ol.format.MVT(),
url: vectorUrl
}),
style: vectorStyle
});
var baseLayer = new ol.layer.Tile({
title: 'osm',
source: new ol.source.XYZ({
url: "http://osm.oslandia.io/styles/osm-bright/{z}/{x}/{y}.png"
})
});
var map = new ol.Map({
target: 'map',
view: new ol.View({
center: [260369,6249082],
zoom: 16
}),
layers: [baseLayer, vectorLayer]
});
var layerSwitcher = new ol.control.LayerSwitcher({
tipLabel: 'Legend', // Optional label for button
groupSelectStyle: 'children' // Can be 'children' [default], 'group' or 'none'
});
map.addControl(layerSwitcher);
var popup = new Popup();
map.addOverlay(popup);
map.on('singleclick', function(evt) {
var info_text = '';
var features = map.getFeaturesAtPixel(evt.pixel);
if (features.length != 0) {
var properties = features[0].getProperties();
info_text = JSON.stringify(properties, null, 2);
}
popup.show(evt.coordinate, '<div><h3>Attributes</h3><p>' + info_text + '</p></div>');
});
</script>
</body>
</html>
```
You can now :
- Try popup
- Change the list of attributes retrieved from the vector tile
- To go further : changer the CSS style of the popup
example of a nicer popup
```javascript
map.on('singleclick', function(evt) {
let infoText = '';
var features = map.getFeaturesAtPixel(evt.pixel);
let i = 0;
for (const feat of features) {
const properties = feat.getProperties();
let featText = '';
for(const prop in properties) {
featText += `<li>${prop}: ${properties[prop]}</li>`;
}
infoText += `
<h4>Feature n°${i}</h4>
<ul>
${featText}
</ul>
`;
i++;
}
popup.show(evt.coordinate, '<div><h3>Attributes</h3>' + infoText + '</div>');
});
```
## Part 5 : PostGIS and pg_tileserv
pg_tileserv exposes layer from tables, but also :
- views
- functions, with parameters
### A vector tile layer from a PostGIS view
We create the following view in PostGIS :
```sql
CREATE OR REPLACE VIEW public.communes_points AS
SELECT
row_number() OVER () AS id
, st_setsrid(c.geom, 2154)::Geometry(Point, 2154) as geom
FROM (
SELECT
(st_dumppoints(
st_generatepoints(
communes.wkb_geometry,
200
, 1)
)
).geom AS geom
, communes.id
FROM
communes
) as c
;
COMMENT ON VIEW public.communes_points IS 'A table of 200 random points per commune (aka arrondissement) in Paris.';
```
We can now :
- Open this layer with `pg_tileserv` in the preview
- Modify parameters ( number of points, seed )
- Add this layer to our application
- Create a new layer as a PostGIS view, displaying buffers around random points in each *arrondissement*
### Another view : urban blocks
From the cadastral parcels, we determine urban blocks : all groups of parcels having common boundaries.
```sql
CREATE view parcelles_blocks AS
SELECT
ogc_fid
, commune::integer * 10000 + ST_ClusterDBSCAN(wkb_geometry, eps := 2, minpoints := 2) over (partition by commune) AS block_id
, wkb_geometry as geom
FROM
parcelles;
```
We can visualize these blocks in QGIS with a **categorized** style.
The parcel table is large, and we can see that the query is quite slow.
We can also aggregate block geometries with the following query :
```sql
SELECT
block_id as id
, st_multi(ST_union(geom))::Geometry(MultiPolygon, 2154) AS geom
, array_agg(ogc_fid) AS parcelles_id
FROM (
SELECT
ogc_fid
, ST_ClusterDBSCAN(wkb_geometry, eps := 2, minpoints := 2) over () AS block_id
, wkb_geometry as geom
FROM
parcelles
) sq
GROUP BY
block_id;
```
We can visualize the result, creating a new table :
```sql
create table blocks as
SELECT
...
;
-- create an index on geometry
CREATE index idx_blocks_geom on blocks using gist(geom);
```
We could also create a view :
```sql
create view blocks_v as
SELECT
...
```
Issues :
- The query is slow
- Impossible to use spatial indexing to filter on the current map viewport, because of the GROUP BY
What we can do :
- Pre-process the data, and use a pre-computed data table instead of a view
- Use a materialized view
- Use a mechanism allowing us to trigger spatial indexing ( see below )
- Change the application logic so that we only work on a data subset, like displaying a single *arrondissement* ( commune ), or around a specific area the user clicked
A solution allowing indexation, left as an exercice for the end of the workshop :
- Create a function using the parameters of the vector tile being generated to introduce in the query a filter evaluating the agregate only on the parcels visible in the viewport
We will see the concept of vector tiles functions next.
### A vector tile layer with filtered data
We will now create a vector tile layer of our parcels, with a filter on a specific Paris *arrondissement*.
We connect to the database and create the following function.
```sql
CREATE OR REPLACE
FUNCTION public.parcelles_filter(
z integer, x integer, y integer,
arrondissement text default '01')
RETURNS bytea
AS $$
WITH
bounds AS (
SELECT ST_TileEnvelope(z, x, y) AS geom
),
mvtgeom AS (
SELECT
ST_AsMVTGeom(ST_Transform(p.wkb_geometry, 3857), bounds.geom) AS geom,
p.id
FROM
parcelles p
, bounds
WHERE
ST_Intersects(p.wkb_geometry, ST_Transform(bounds.geom, 2154))
AND
p.commune LIKE ('%' || arrondissement)
)
SELECT
ST_AsMVT(mvtgeom, 'public.parcelles_filter')
FROM
mvtgeom;
$$
LANGUAGE 'sql'
STABLE
PARALLEL SAFE;
COMMENT ON FUNCTION public.parcelles_filter IS 'Filter the parcel table with the parameter "arrondissement" on the commune field.';
```
We can now test pg_tileserv :
- Connect on pg_tileserv instance
- visualize the new data layer
- input a parameter to visualize the result
### A vector tile layer with processed data
We will now create a data layer consisting in the result of a geometric processing on the data, with multiple parameters : the intersection of our parcels for a buffor of `dist` around a `lat`, `lon` point.
```sql
CREATE OR REPLACE
FUNCTION public.parcelles_radius(
z integer, x integer, y integer,
click_lon float8 default 2.35,
click_lat float8 default 48.86,
radius float8 default 100)
RETURNS bytea
AS $$
WITH
args AS (
SELECT
ST_TileEnvelope(z, x, y) AS bounds,
ST_Transform(ST_SetSRID(ST_MakePoint(click_lon, click_lat), 4326), 2154) AS click
),
mvtgeom AS (
SELECT
ST_AsMVTGeom(
ST_Transform(
ST_Intersection(
p.wkb_geometry,
ST_Buffer(args.click, radius)),
3857),
args.bounds) AS geom,
p.id
FROM
parcelles p
, args
WHERE
ST_Intersects(p.wkb_geometry, ST_Transform(args.bounds, 2154))
AND
ST_DWithin(p.wkb_geometry, args.click, radius)
LIMIT 10000
)
SELECT ST_AsMVT(mvtgeom, 'public.parcelles_radius') FROM mvtgeom
$$
LANGUAGE 'sql'
STABLE
PARALLEL SAFE;
COMMENT ON FUNCTION public.parcelles_radius IS E'From a given point of coordinates ( click_lon, click_lat ) and a given radius, return all parcels intersecting this circle, with a geometry cut at the corresponding radius.';
```
We can now verify this new layr and test it with `pg_tileserv`.
- Open pg_tileserv homepage
- Test with various values
### Integrate into our web application
We can now integrate this function into our web application.
Edit `index4.html` with the following content :
```htmlmixed=
<!doctype html>
<html lang="en">
<head>
<meta charset="utf-8">
<meta name="viewport" content="width=device-width, initial-scale=1">
<title>OpenLayers Vector Tiles</title>
<!-- CSS/JS for OpenLayers map -->
<link rel="stylesheet" href="https://openlayers.org/en/v6.1.1/css/ol.css" type="text/css">
<script src="https://openlayers.org/en/v6.1.1/build/ol.js"></script>
<script src="https://unpkg.com/ol-layerswitcher@3.7.0"></script>
<script src="https://unpkg.com/ol-popup@4.0.0"></script>
<style>
body {
padding: 0;
margin: 0;
}
html, body, #map {
height: 100%;
width: 100%;
font-family: sans-serif;
}
#meta {
background-color: rgba(255,255,255,0.75);
color: black;
z-index: 2;
position: absolute;
top: 10px;
left: 20px;
padding: 10px 20px;
margin: 0;
}
</style>
<link rel="stylesheet" href="https://unpkg.com/ol-layerswitcher@3.7.0/src/ol-layerswitcher.css" />
<link rel="stylesheet" href="https://unpkg.com/ol-popup@4.0.0/src/ol-popup.css" />
</head>
<body>
<div id="meta">
<h2>Our cartographic application !</h2>
<ul>
<li><a href="https://localhost:7800">pg_tileserv server</a></li>
<li>
<form>
<label>Radius :</label>
<input type="number" id="radius" name="radius" min="10" max="2000" value="100" step="50">
</form>
</li>
</ul>
</div>
<div id="map"></div>
<script>
var vectorServer = "http://localhost:7800/";
var vectorSourceLayer = "public.parcelles";
var vectorProps = "?properties=ogc_fid,id,commune,prefixe,section,numero,contenance,created,updated"
var vectorUrl = vectorServer + vectorSourceLayer + "/{z}/{x}/{y}.pbf" + vectorProps;
var vectorStyle = new ol.style.Style({
stroke: new ol.style.Stroke({
width: 2,
color: "#ff00ff99"
}),
fill: new ol.style.Fill({
color: "#ff00ff33"
})
});
var vectorStyleBlue = new ol.style.Style({
stroke: new ol.style.Stroke({
width: 1,
color: "#0000ff99"
}),
fill: new ol.style.Fill({
color: "#0000ff33"
})
});
var vectorLayer = new ol.layer.VectorTile({
title:'Parcelles'
,source: new ol.source.VectorTile({
format: new ol.format.MVT(),
url: vectorUrl
}),
style: vectorStyle
});
var parcelles_radius = new ol.layer.VectorTile({
title:'Parcelles radius'
,source: new ol.source.VectorTile({
format: new ol.format.MVT(),
url: vectorServer + "public.parcelles_radius" + "/{z}/{x}/{y}.pbf" + "?properties=id"
}),
style: vectorStyleBlue
});
var baseLayer = new ol.layer.Tile({
title: 'osm',
source: new ol.source.XYZ({
url: "http://osm.oslandia.io/styles/osm-bright/{z}/{x}/{y}.png"
})
});
var map = new ol.Map({
target: 'map',
view: new ol.View({
center: [260369,6249082],
zoom: 16
}),
layers: [baseLayer, vectorLayer, parcelles_radius]
});
var layerSwitcher = new ol.control.LayerSwitcher({
tipLabel: 'Legend', // Optional label for button
groupSelectStyle: 'children' // Can be 'children' [default], 'group' or 'none'
});
map.addControl(layerSwitcher);
var popup = new Popup();
map.addOverlay(popup);
map.on('singleclick', function(evt) {
var click = ol.proj.transform(evt.coordinate, 'EPSG:3857', 'EPSG:4326');
var radius = document.getElementById('radius').value;
var source = parcelles_radius.getSource();
source.setUrl(vectorServer
+ "public.parcelles_radius"
+ "/{z}/{x}/{y}.pbf"
+ "?properties=id"
+ "&click_lon=" + click[0]
+ "&click_lat=" + click[1]
+ "&radius=" + radius)
source.refresh({force: true});
});
</script>
</body>
</html>
```
## Part 6 : PostGIS and pg_featureserv
### new parcelles_radius
The last function `public.parcelle_radius` might be better served by pg_featureserv, as it gives back a set of features.
- create a `postgisftw` schema to host our function
- adapt the previous function to be suitable for pg_featureserv (see official documentation).
```sql
CREATE OR REPLACE FUNCTION postgisftw.parcelles_radius(
click_lon double precision DEFAULT 2.35,
click_lat double precision DEFAULT 48.86,
radius double precision DEFAULT 100)
RETURNS TABLE(wkb_geometry geometry, id character varying, numero character varying, contenance numeric)
LANGUAGE sql
STABLE PARALLEL SAFE
AS $function$
WITH
args AS (
SELECT
ST_Transform(ST_SetSRID(ST_MakePoint(click_lon, click_lat), 4326), 2154) AS click
)
SELECT
ST_Transform(
ST_Intersection(
p.wkb_geometry,
ST_Buffer(args.click, radius)),
3857) as wkb_geometry,
p.id, p.numero, p.contenance
FROM
parcelles p
, args
WHERE
ST_DWithin(p.wkb_geometry, args.click, radius)
LIMIT 10000
$function$;
COMMENT ON FUNCTION postgisftw.parcelles_radius IS E'From a given point of coordinates ( click_lon, click_lat ) and a given radius, return all parcels intersecting this circle, with a geometry cut at the corresponding radius.';
```
- check the home page of pg_featureserv to see if the function has been loaded correctly
- edit our index4.html to use this function instead of pg_tileserv url
```htmlmixed=
<!doctype html>
<html lang="en">
<head>
<meta charset="utf-8">
<meta name="viewport" content="width=device-width, initial-scale=1">
<title>OpenLayers Vector Tiles</title>
<!-- CSS/JS for OpenLayers map -->
<link rel="stylesheet" href="https://openlayers.org/en/v6.1.1/css/ol.css" type="text/css">
<script src="https://openlayers.org/en/v6.1.1/build/ol.js"></script>
<script src="https://unpkg.com/ol-layerswitcher@3.7.0"></script>
<script src="https://unpkg.com/ol-popup@4.0.0"></script>
<style>
body {
padding: 0;
margin: 0;
}
html, body, #map {
height: 100%;
width: 100%;
font-family: sans-serif;
}
#meta {
background-color: rgba(255,255,255,0.75);
color: black;
z-index: 2;
position: absolute;
top: 10px;
left: 20px;
padding: 10px 20px;
margin: 0;
}
</style>
<link rel="stylesheet" href="https://unpkg.com/ol-layerswitcher@3.7.0/src/ol-layerswitcher.css" />
<link rel="stylesheet" href="https://unpkg.com/ol-popup@4.0.0/src/ol-popup.css" />
</head>
<body>
<div id="meta">
<h2>Our cartographic application !</h2>
<ul>
<li><a href="https://localhost:7800">pg_tileserv server</a></li>
<li>
<form>
<label>Radius :</label>
<input type="number" id="radius" name="radius" min="10" max="2000" value="100" step="50">
</form>
</li>
</ul>
</div>
<div id="map"></div>
<script>
var vectorServer = "http://localhost:7800/";
var vectorSourceLayer = "public.parcelles";
var vectorProps = "?properties=ogc_fid,id,commune,prefixe,section,numero,contenance,created,updated"
var vectorUrl = vectorServer + vectorSourceLayer + "/{z}/{x}/{y}.pbf" + vectorProps;
var vectorStyle = new ol.style.Style({
stroke: new ol.style.Stroke({
width: 2,
color: "#ff00ff99"
}),
fill: new ol.style.Fill({
color: "#ff00ff33"
})
});
var vectorStyleBlue = new ol.style.Style({
stroke: new ol.style.Stroke({
width: 1,
color: "#0000ff99"
}),
fill: new ol.style.Fill({
color: "#0000ff33"
})
});
var vectorLayer = new ol.layer.VectorTile({
title:'Parcelles'
,source: new ol.source.VectorTile({
format: new ol.format.MVT(),
url: vectorUrl
}),
style: vectorStyle
});
var parcelles_radius = new ol.layer.Vector({
title:'Parcelles radius'
,source: new ol.source.Vector({
format: new ol.format.GeoJSON(),
url: "http://localhost:9000/functions/parcelles_radius/items.json"
}),
style: vectorStyleBlue
});
var baseLayer = new ol.layer.Tile({
title: 'osm',
source: new ol.source.XYZ({
url: "http://osm.oslandia.io/styles/osm-bright/{z}/{x}/{y}.png"
})
});
var map = new ol.Map({
target: 'map',
view: new ol.View({
center: [260369,6249082],
zoom: 16
}),
layers: [baseLayer, vectorLayer, parcelles_radius]
});
var layerSwitcher = new ol.control.LayerSwitcher({
tipLabel: 'Legend', // Optional label for button
groupSelectStyle: 'children' // Can be 'children' [default], 'group' or 'none'
});
map.addControl(layerSwitcher);
var popup = new Popup();
map.addOverlay(popup);
map.on('singleclick', function(evt) {
var click = ol.proj.transform(evt.coordinate, 'EPSG:3857', 'EPSG:4326');
var radius = document.getElementById('radius').value;
var source = parcelles_radius.getSource();
source.setUrl("http://localhost:9000/functions/parcelles_radius/items.json"
+ "?properties=id"
+ "&click_lon=" + click[0]
+ "&click_lat=" + click[1]
+ "&radius=" + radius
+ "&limit=" + 2000)
source.refresh({force: true});
});
</script>
</body>
</html>
```
<details>
<summary><strong>Correction</strong> Steps to follow</summary>
<ul>
<li>change the `parcelles_radius` layer type to Vector, with a vector source and a Geojson format</li>
<li>change the click handler to use the pg_featureserv function instead of the pg_tileserv url</li>
</ul>
</details>
- What do you observe?
:::info
<details><summary>Answer</summary>
- there is only one web request done
- and it is much faster
- why?
- pg_tileserv need to query all the visible tiles, even though most of them would contain no data, and it does one query for each of these tiles
- it needs to do more transformation
- there are 2 geometric filters, but only one of them can use the index
</details>
:::
- What are your conclusion on pg_featureserv and pg_tileserv?
## More exercices and ideas
We now have all elements in place to add new features, coupling the web Javascript interface directly with the cartographic server pg_tileserv, leveraging PostGIS capabilities.
A few feature ideas :
- A parcel search feature, given a parcel ID
- a building search feature, given a name
- A parcel layer including a computed total building surface as properties
- with a specific style for parcels with a built surface over 80%
- A layer of parcels with building footprints removed from the geometry
- A voronoï polygon layer of parcels
- Optimizing the block creation queries and layers
# Links and resources
Some linkes and resources :
- [pg_tileserv documentation](https://access.crunchydata.com/documentation/pg_tileserv/latest/)
- [OpenLayers documentation](https://openlayers.org/en/latest/doc/)
- [Layer Switcher](https://github.com/walkermatt/ol-layerswitcher)
- [Popup](https://github.com/walkermatt/ol-popup)
- [Workshop OpenLayers](https://openlayers.org/workshop/en/)
- [OL extensions](https://viglino.github.io/ol-ext/)
- [PostGIS Documentation - référence chapter](https://postgis.net/docs/reference.html)
# Notes on the virtual machine
The virtual Machine is a Virtualbox instance with Ubuntu 20.04 installed, and the following packages added :
```bash=
sudo -s
echo "deb https://qgis.org/ubuntu-nightly-release groovy main" > /etc/apt/sources.list.d/qgis.list
wget -4 -qO - https://qgis.org/downloads/qgis-2020.gpg.key | sudo gpg --no-default-keyring --keyring gnupg-ring:/etc/apt/trusted.gpg.d/qgis-archive.gpg --import
chmod a+r /etc/apt/trusted.gpg.d/qgis-archive.gpg
apt update
apt install postgis qgis postgresql-12 postgresql-12-postgis-3 postgresql-12-postgis-3-scripts
```
To configure postgres password :
```bash=
sudo -u postgres psql
# alter user postgres with password 'oslandia';
```