Introduction

In my previous post I looked at how to turn a Lidar point cloud into a simple canopy height model using freely available tools. https://rpubs.com/dgolicher/lidar

A problem with raw lidar data is that there is simply too much of it at too fine a resolution for regional analysis. The raw data can be held on servers as tiles and processed by powerful servers. but even so fine scale raster data processing quickly becomes complex. The conventional solution is to derive statistics based on the fine grained rasters onto coarser resolution grids. However calculating a single statistic such as mean canopy height loses too much information. It is usually necessary to retain more information on the distribution of the fine grained data.

One advantage of turning raster layers into vector representations is that they can be easily filtered and thus reduced in size still more. It can be especially useful to add these layers to spatial data bases.

In this tutorial I will give a simple example of how data can be processed to a local PostgIS data base. In a production environment the results would be added to a regional data layer.

Setting up the database

The code as given will only work directly on Linux, as it calls system commands from R. In order to set up a local PostGIS system the instructions given here need to be followed first.

https://rpubs.com/dgolicher/6373

I set up a database called “lidar” with the createdb command. From R this works, providing the permissions are all OK.

system("dropdb lidar")
system("createdb lidar")

I then setup an ODBC connection to this data base as shown in the tutorial linked above. My odbc.ini file for my 64 bit Ubuntu laptop looks like this.

[lidar]
Driver = /usr/lib/x86_64-linux-gnu/odbc/psqlodbcw.so
Database = lidar
Servername = localhost
Username = postgres
Password = postgres
Protocol = 8.2.5
ReadOnly = 0

Then from R the extensions can be added.

library(RODBC)
con <- odbcConnect("lidar")
odbcQuery(con, "CREATE EXTENSION POSTGIS")
odbcQuery(con, "CREATE EXTENSION PLR")

Loading raster layer from QGIS

Anyone familiar with PostGIS and QGIS will know that vector layers can be added to the data base with a simple drag and drop. However there is (at time of writing 27/11/2016) still no such simple way of loading a raster layer from the canvas. Raster layers are loaded using a raster2psql command.

By pasting some of the paramters into the command (file name and grid size) in R a command can be put together to run on the system. The script below can be added to the QGIS toolbox as an R script raster loader to the (hardwired) lidar data base. It is simple to change it to include a prarmeter for choosing another data base.

##Raster processing=group
##Raster = raster
##grid_side= number 10
library(raster)
library(RODBC)
r<-Raster
r_nm<-r@file@name
command <- paste("raster2pgsql -d  -M -R ",r_nm, " -F -t ,grid_side,"x",grid_side," tmp|psql -d lidar",sep="")
system(command)

The script produces the following dialogue in QGIS.

Running it loads the CHM layer into the database as an external (out of deatabase) raster layer.

Aggregating the values to a grid

The first step has already set up a potential grid to use, as the data are tiled. There are functions in PostGIS to output the geometry of the tiles as polygons and the values. All that is needed is to setup some PLR functions to calculate the stats we need and write a SQL query to form a new table.

So to add a function to calulate the median to the lidar database from R we first connect to it using ODBC then run the query.

library(RODBC)
con<-odbcConnect("lidar")
query<-"CREATE OR REPLACE FUNCTION median (float[]) RETURNS float AS '
x<-arg1
x<-as.numeric(as.character(x))
x<-na.omit(x)
median(x)'
LANGUAGE 'plr' STRICT;"
odbcQuery(con,query)

Other functions for calculating statistics can follow the same pattern. So this query adds functions for quantiles to the median.

query<-"CREATE OR REPLACE FUNCTION median (float[]) RETURNS float AS '
x<-arg1
x<-as.numeric(as.character(x))
x<-na.omit(x)
median(x)'
LANGUAGE 'plr' STRICT;
CREATE OR REPLACE FUNCTION q10 (float[]) RETURNS float AS '
x<-arg1
x<-as.numeric(as.character(x))
x<-na.omit(x)
quantile(x,0.1,na.rm=TRUE)'
LANGUAGE 'plr' STRICT;
CREATE OR REPLACE FUNCTION q90 (float[]) RETURNS float AS '
x<-arg1
x<-as.numeric(as.character(x))
x<-na.omit(x)
quantile(x,0.9,na.rm=TRUE)'
LANGUAGE 'plr' STRICT;"

odbcQuery(con,query)

Now all we need is to execute a query to produce these aggregated statistics and add them to a vector grid based on the chosen tile dimensions.

query<-"drop table if exists tmp2;
create table tmp2 as
select rid, st_envelope(rast) geom,
q10((st_dumpvalues(rast)).valarray) q10,
median((st_dumpvalues(rast)).valarray) median,
q90((st_dumpvalues(rast)).valarray) q90
from tmp"

odbcQuery(con,query)

The full script

Adding the following code to the QGIS R toolbox should allow any raster layer to be loaded to the lidar data base as a table called “tmp” and then aggregated using the three functions. In general PostGIS runs faster than the equvalent code in native R, although there may be some overhead in loading to the database first. One advantage for desktop use is that there is no need to produce a graticule layer first. In a production seting the aim would be to combine the results to form a queryable layer at a regional scale.

##Raster processing=group
##Raster = raster
##grid_side= number 10
library(raster)
library(rgdal)
library(RODBC)
r<-Raster
r_nm<-r@file@name

command <- paste("raster2pgsql -d  -M -R ",r_nm, " -F -t ",grid_side,"x",grid_side," tmp|psql -d lidar",sep="")
system(command)

con<-odbcConnect("lidar")

query<-"CREATE OR REPLACE FUNCTION median (float[]) RETURNS float AS '
x<-arg1
x<-as.numeric(as.character(x))
x<-na.omit(x)
median(x)'
LANGUAGE 'plr' STRICT;
CREATE OR REPLACE FUNCTION q10 (float[]) RETURNS float AS '
x<-arg1
x<-as.numeric(as.character(x))
x<-na.omit(x)
quantile(x,0.1,na.rm=TRUE)'
LANGUAGE 'plr' STRICT;
CREATE OR REPLACE FUNCTION q90 (float[]) RETURNS float AS '
x<-arg1
x<-as.numeric(as.character(x))
x<-na.omit(x)
quantile(x,0.9,na.rm=TRUE)'
LANGUAGE 'plr' STRICT;"

odbcQuery(con,query)

query<-"drop table if exists tmp2;
create table tmp2 as
select rid, st_envelope(rast) geom,
q10((st_dumpvalues(rast)).valarray) q10,
median((st_dumpvalues(rast)).valarray) median,
q90((st_dumpvalues(rast)).valarray) q90
from tmp"
odbcQuery(con,query)

The script can be downloaded from here.

https://dgolicher.bitbucket.io/lidar/qgis_aggregate_script.rsx

Layer added to canvas from PostGIS