10.2 Geospatial data processing in SQL - Part 1¶

Tip

When executing SQL statements in QGIS DB Manager check also Log messages window (PostGIS tab) for possible SQL error messages.

Metadata tables¶

Geometry_columns¶

Table geometry_columns describes the available feature tables and their geometry properties.

select * from geometry_columns;
select type,count(*) from geometry_columns group by type;
select distinct srid from geometry_columns;

Spatial_ref_sys¶

Table spatial_ref_sys contains coordinate reference system (CRS) definitions.

select * from spatial_ref_sys;
select * from spatial_ref_sys where srid = 4326;

Let's combine information from geometry_columns and spatial_ref_sys:

select g.f_table_name,g.type,s.srtext from geometry_columns as g
join spatial_ref_sys as s
using (srid);

Geometry vs Geography¶

Geometry area is computed in map units:

select sum(st_area(geom)) from landuse_a where fclass = 'vineyard';
select sum(st_area(st_transform(geom, 3035)))/1e4 from landuse_a where fclass = 'vineyard';

Geography represents feature in geodetic coordinate systems which models the earth using an ellipsoid. Geography calculations are in meters.

select sum(st_area(geom::geography))/1e4 from landuse_a where fclass = 'vineyard';

Note

Compare computed area on an ellipsoid and projected CRS (EPSG 3035 and 5514):

select sum(st_area(geom::geography))/1e4 as area from landuse_a where fclass = 'vineyard'
union all
select sum(st_area(st_transform(geom, 3035)))/1e4 from landuse_a where fclass = 'vineyard'
union all
select sum(st_area(st_transform(geom, 5514)))/1e4 from landuse_a where fclass = 'vineyard';

Spatial operators¶

BBox operators¶

This spatial operators operate on bounding boxes:

&& — Returns TRUE if A's 2D bounding box intersects B's 2D bounding box.
@ — Returns TRUE if A's bounding box is contained by B's.
~= — Returns TRUE if A's bounding box is the same as B's.
...

Full list: https://postgis.net/docs/reference.html#operators-bbox

select
st_geomfromtext('POLYGON((1 0, 2 1, 1 2, 0 1, 1 0))')
&&
st_geomfromtext('POLYGON((3 2, 4 3, 3 4, 2 3, 3 2))');

Task

Try different bbox operators: &<, &>, <<, ...

Now let's focus on OSM data (vineyards):

create view vineyards as select geom from landuse_a where fclass = 'vineyard';

and apply && spatial operator on roads:

select r.* from roads as r
join vineyards as v
on r.geom && v.geom;

And display result in QGIS.

Tip

Take a look on query planner (explain select ... and explain analyze select ...). To improve performence create index on landuse_a.fclass.

create index on landuse_a (fclass);

Distance operators¶

<-> — Returns the 2D distance between A and B.
<#> — Returns the 2D distance between A and B bounding boxes.
...

Full list: https://postgis.net/docs/reference.html#operators-distance

select
st_geomfromtext('POLYGON((1 0, 2 1, 1 2, 0 1, 1 0))')
<->
st_geomfromtext('POLYGON((3 2, 4 3, 3 4, 2 3, 3 2))');

Task

Compare with <#> distance operator.

Let's find the nearest library (current position: https://mapy.cz/s/coremubako):

with cp as 
(select st_setsrid(st_point(14.3881100, 50.1041200), 4326)::geography as geom) 
select name,p.geom::geography <-> cp.geom as dist from pois as p
cross join cp
where fclass = 'library'
order by dist
limit 5;

Note

Compare with st_distance().

Spatial filters¶

Spatial relationships have been explain in the previous lesson.

Examples above:

select st_intersects(
st_geomfromtext('POLYGON((1 0, 2 1, 1 2, 0 1, 1 0))'),
st_geomfromtext('POLYGON((3 2, 4 3, 3 4, 2 3, 3 2))'));

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

To get answer how spatial filtering is performed, examine query planner by explain select ....

Note

Query planner gives also explanation why distinct is needed here.

Spatial filtering is performed in two steps:

(1) Select geometries using bbox operator (&&);

(2) Apply spatial relationship predicate only on subset of geometries from the step above.

Task

Try various spatial predicates (st_within(), st_contains(), st_touches()).

Example:

select distinct ww.* from waterways as ww
join water_a as w
on st_touches(ww.geom, w.geom);

Let's take a closer look on spatial filter applied in lesson 9: find nearest cinema based on our current location. This task may be elegantly (and more precisely - buffer is approximated by a polygon - note that rings are represented by rings/linestrings not curves) solved by st_dwithin().

with cinema as
(select st_transform(geom, 3035) as geom, name from pois where fclass = 'cinema'),
aoi as
(select st_transform(st_setsrid(st_point(14.3886592, 50.1042994), 4326), 3035) as geom)
select c.* from cinema as c
join aoi
on st_dwithin(c.geom, aoi.geom, 3000);

Task

Compare result with

with cinema as
(select st_transform(geom, 3035) as geom, name from pois where fclass = 'cinema'),
aoi as
(select st_buffer(st_transform(st_setsrid(st_point(14.3886592, 50.1042994), 4326), 3035), 3000) as geom)
select c.* from cinema as c
join aoi
on st_within(c.geom, aoi.geom);

and explain difference. Try to increase buffer precision, see PostGIS documentation.