10.3 Geospatial data processing in SQL - Part 2¶

LUCAS Data¶

In this lesson we will combine OpenStreetMap data with LUCAS dataset. LUCAS (Land Use and Coverage Area frame Survey, https://ec.europa.eu/eurostat/web/lucas) is an activity that has performed repeated in situ surveys over Europe every three years since 2006. Dataset identifies changes in land use and land cover in the EU.

Let's download LUCAS data (year 2018) for selected country in CSV format: https://ec.europa.eu/eurostat/web/lucas/data/primary-data/2018

CSV file can be imported into PostGIS database by DB Manager Import Layer/File:

Note

Don't forget to convert field names to lowercase

Create geometry¶

Geometry constructors have been introduced in previous lesson.

Add a new geometry column:

alter table lucas_cz_2018 add column geom geometry(point, 4326);

and create geometry by st_point() function:

update lucas_cz_2018 set geom=st_setsrid(st_point(cast(gps_long as float), cast(gps_lat as float)), 4326);

Check result by SQL:

select count(*) from lucas_cz_2018 where st_isvalid(geom) = false;
select count(*) from lucas_cz_2018 where st_isempty(geom);
select st_x(geom), st_y(geom) from lucas_cz_2018;

Drop geometry for non-sense locations:

update lucas_cz_2018 set geom = NULL where gps_prec = '8888';

Add spatial index:

create index on lucas_cz_2018 using gist (geom);

Load result of SQL query into map canvas:

Compare land use in LUCAS and OSM¶

Select LUCAS points located in vineyards:

with vineyard as
(select geom from landuse_a where fclass = 'vineyard')
select l.point_id,l.lu1 from lucas_cz_2018 as l 
join vineyard as v
on st_intersects(l.geom, v.geom)
where l.geom is not NULL;

Compare with LUCAS LU1 attribute with landuse OSM attribute:

select l.point_id,l.lu1,o.fclass from lucas_cz_2018 as l 
join landuse_a as o
on st_intersects(l.geom, o.geom)
where l.geom is not NULL;

Note

You may import LUCAS LU1 coding table in CSV format into PostGIS and get comparision of class names.

with lu as
(
 select l.point_id,l.lu1,o.fclass from lucas_cz_2018 as l 
 join landuse_a as o
 on st_intersects(l.geom, o.geom)
 where l.geom is not NULL)
select lu.point_id,lu1.name as lu_lucas,lu.fclass as lu_osm from lu
join lu1
on lu.lu1 = lu1.code;

Spatial overlay¶

Let's continue in our example with vineyards:

select st_intersection(r.geom, v.geom) from roads as r
join vineyards as v
on st_intersects(r.geom, v.geom);

Example from previous lesson can be significantly simplified in PostGIS:

create table natural_areas as
select l.osm_id, st_intersection(nuts.geom, st_transform(l.geom, 3035)) intersected_geom
from landuse_a as l, nuts
where fclass = 'nature_reserve'
and nuts.nuts_id = 'CZ010'
and st_intersects(nuts.geom, st_transform(l.geom, 3035)) = True;

Note

You have to import nuts data into PostGIS database. The simpliest way is to use DB Manager (Import Layer/File).

Primary key and spatial index may be created by alter table and create index statements:

alter table natural_area add primary key(osm_id);

create index on natural_area using gist (intersected_geom);

Let's take a closer look at geometries created by intersection of landuse and nuts features:

select distinct st_geometrytype(intersected_geom) from natural_areas;
select st_geometrytype(intersected_geom) as gtype, count(*) from natural_areas group by gtype;
select st_numgeometries(intersected_geom) from natural_areas where st_geometrytype(intersected_geom) = 'ST_MultiPolygon';

Mixture of muplipolygons and polygons may be not desired. Let's split multipolygon features into multiple polygon features by st_dump():

select (st_dump(intersected_geom)).geom from natural_areas;