Database Setup

How to install, setup and populate the spatial database

1. Install the latest version of postgresql with the postgis extension for your system then install pgpointcloud

apt install postgresql-10-postgis-2.4 postgis postgresql-server-dev-10

2. Create an user for the database and a db assosciated to id:

First log in as the postgres root user

sudo su - postgres

Then create a new role and grant it the superuser privileges (needed to create extensions):

createuser --pwprompt -s dboperator

Then create the db:

createdb -O dboperator terrain_ans

3. Install pgpointcloud by following the project wiki: https://github.com/pgpointcloud/pointcloud
4. Create the lidar tables and the extensions in the postgis db:

CREATE EXTENSION postgis;
CREATE EXTENSION pointcloud;
CREATE EXTENSION pointcloud_postgis;
CREATE TABLE lidar_toscana (
    id SERIAL PRIMARY KEY,
    pa PCPATCH(1)
);
CREATE INDEX lidar_toscana_index ON lidar_toscana USING GIST(PC_EnvelopeGeometry(pa));
CREATE TABLE lidar_lyon (
    id SERIAL PRIMARY KEY,
    pa PCPATCH(1)
);

4. Install pdal and it's dependencies (laszip support mandatory for certain dataset)

apt install pdal

5. Download the lidar files for your zone and add them using the following script

#!/bin/bash

for file in `ls  ~/data/Toscana/lidar/asc/*.asc`; do
	pdal pipeline toscana.json --readers.gdal.filename=$file && mv $file $file.loaded
	echo $file loaded
done 

The file toscana.json contains the example pipeline for the Tuscanian lidar data. It must be changed for different sources

{
  "pipeline":[
    {
      "type" : "readers.gdal"
    },
    {
     "type":"filters.reprojection",
     "in_srs":"EPSG:3003",
     "out_srs":"EPSG:4326"
    },
    {
     "type" : "filters.ferry",
     "dimensions": "band-1=Z"
    },
    {
     "type":"filters.chipper",
     "capacity":"600"
    },
    {
     "type": "writers.pgpointcloud",
     "connection":"host='hostname' dbname='db_name' user='user' password='password'",
     "table" : "lidar_toscana",
     "column": "pa",
     "srid": "4326",
     "compression":"dimensional",
     "output_dims": "X,Y,Z"
    }
  ]
}

6. Download the osm building files from an osm mirror
7. Insert the osm building using the shp2pgsql command (similar flow can be used for other vectorial data)

shp2pgsql -s 4326 -I ~/osm/lyon/gis.osm_buildings_a_free_1.shp  lyon_buildings | sudo -u postgres psql
shp2pgsql -s 4326 -I ~/osm/centro_italy/gis.osm_buildings_a_free_1.shp  centro_buildings | sudo -u postgres psql

Create the following helper postgres functions:

CREATE OR REPLACE FUNCTION public.st_buffer_meters(
    geometry,
    double precision)
  RETURNS geometry AS
$BODY$
DECLARE
orig_srid int;
utm_srid int;
 
BEGIN
orig_srid:= ST_SRID($1);
utm_srid:= utmzone(ST_Centroid($1));
 
RETURN ST_transform(ST_Buffer(ST_transform($1, utm_srid), $2), orig_srid);
END;
$BODY$
  LANGUAGE plpgsql IMMUTABLE
  COST 100;

and

CREATE OR REPLACE FUNCTION public.utmzone(geometry)
  RETURNS integer AS
$BODY$
DECLARE
geomgeog geometry;
zone int;
pref int;
 
BEGIN
geomgeog:= ST_Transform($1,4326);
 
IF (ST_Y(geomgeog))>0 THEN
pref:=32600;
ELSE
pref:=32700;
END IF;
 
zone:=floor((ST_X(geomgeog)+180)/6)+1;
 
RETURN zone+pref;
END;
$BODY$
  LANGUAGE plpgsql IMMUTABLE
  COST 100;