Building utility functions in PostGIS

We want the data that we place in our PostGIS database to be easily extracted by students and technical staff with limited knowledge of spatial SQL. Some queries, particularly those involving raster data can be very daunting. We could try to use a mapserver or build queries for each user. However we would like to keep things as generic and flexible as possible.

Most of the common queries involve localities, climate and land use. So it is worth building some simple extractor functions that take the geometry of any layer as an argument and return the information that is most useful.

Here is an example based on a post over a year ago that finds the nearest locality to any point in a layer which has a population over a certain cut off. Things have changed a bit since PostGIS1.5 and we can now use geography based functions to calculate distances and areas.

A simple trick is to first build a view as a template for a function. This not only allows you to test out the function first but if you don’t drop the view you are left with a recordset that has registered the variables in the right form that can be returned from the function.

We want to find the closest green points to the red points given some lower threshold of population.

So we can test a query first as a view.


create view mas_cerc as
select distinct on (d.geom) d.nom_loc,dist_km from
(select h.id, nom_loc, h.geom,
round(cast(st_distance_sphere(w.geom,h.geom)/1000 as numeric),3) as
dist_km from chiapas.gbif h,
(select * from chiapas.adm_loc2010 where pobtot>100) w
where st_dwithin(w.geom,h.geom,0.2)
order by h.geom,dist_km) as d

This is more complex than a student would be expected to build. The view uses two tables. One contains the key information that will be extracted (chiapas.adm_loc2010). However we may want to extract this information for any other point table in the future. So we want to build a function that just takes the geometry of a table, but is also takes the population cut off as an argument.


create function loc_cerc(geometry,float) returns setof mas_cerc as
'select distinct on (d.geom) d.nom_loc,dist_km from
(select nom_loc, $1 as geom,
round(cast(st_distance_sphere(w.geom,$1)/1000 as numeric),3) as
dist_km
from
(select * from chiapas.adm_loc2010 where pobtot>$2) w
where st_dwithin($1,w.geom,0.2)
order by $1,dist_km) as d'
language sql;

Now we can get the nearest localities to all the oaks simply by writing


select *,

(loc_cerc(geom,100)).*

from chiapas.gbif where genus like 'Quercus'

Notice that as the function returns a recordset you need to surround it with brackets and add ().* to get a column for each return value.

Check that it does what it should be filtering the localities layer to only those over 100 and labelling the points with contrasting colours.

It works. So now all that is needed are three lines of SQL to use it with any point data set within Chiapas.

A similar situation occurs with the Worldclim data. There are two issues here. The first is that users cannot import the raster layers themselves into their own PostGIS data bases very easily. The second is that Worldclim potentially provides more information that is usuallly needed. So this function will simplify the output.


create view bioclim1 as
select
st_value(rast,1,geom) Elev ,
st_value(rast,2,geom)/10 TempAnnual,
st_value(rast,7,geom)/10 MesMasFrio ,
st_value(rast,6,geom)/10 MesMasCalido ,
st_value(rast,12,geom) PrecAnual

from
(select * from chiapas.gbif where genus like 'Quer%') treestab,
rasters.bioclim10x10 where st_intersects(rast,geom)

Using the same trick to build the function.


create function bioclim(geometry) returns setof bioclim1 as
'select
st_value(rast,1,$1) Elev ,
st_value(rast,2,$1)/10 TemAnnual,
st_value(rast,7,$1)/10 MesMasFio ,
st_value(rast,6,$1)/10 MesMasCalido ,
st_value(rast,12,$1) PrecAnual

from

rasters.bioclim10x10
where st_intersects(rast,$1)'

language sql;

The view and the associated function can of course be easily extended to pull out more of the bioclim variables if they are wanted.

So now all we need to do is write

select *,
(bioclim(geom)).*
from chiapas.gbif where genus like 'Quercus';

To get the result we want.

So anyone logged into the institutional server can use this function. If they contribute their own data to the database they can run these, and more utility functions on them.

A final utility. Producing a list of all the species found around a point.


create view sp_list as
select * from
(select distinct on (g.nombre) family, genus,species
from chiapas.gbif g,
(select * from chiapas.adm_loc2010 where cve_loc like '0001') as l
where st_dwithin(l.geom,g.geom,0.1)) a
order by family,genus

create function spec_list(geometry,float) returns setof sp_list as
'select * from
(select distinct on (g.nombre) family, genus,species
from chiapas.gbif g
where st_dwithin($1,g.geom,$2)) a
order by family,genus'
language sql;

So to get a list of all the species within a 0.1 degree radius of San Cristobal we just need to write.


select nom_loc,
(spec_list(geom,0.1)).*
from chiapas.adm_loc2010 where nom_loc like 'SAN CRISTÓBAL DE LAS CASAS%'

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s