Maths, geo and computer science ...

December 6, 2017

Overview of the problem

We have two layers :

• A layer of the different Provinces of Laos
• A layer of the villages of Laos

The problem is the following: How can I request all the villages that are inside a given Province ?

TL;DR: Using PostGIS and the statement ST_Contains.

Overview of the data

The two layers used can be found on the OpenStreetMap website for laos. You can download the zipfile here. In order to visualize the data, we can use QGIS. Just click on Layer » Add a Vector Layer and select the zip file, then select the .shp files in the next window. You should see the following map.

Let’s focus on the Province of Vientiane and list all the places inside of it.

You can see the table of attributes for the two layers, the places have a column name with names in english and lao.

PostGIS is the spatial extension of PostgreSQL. If you have not install it, instructions are here.

Once installed, let’s create a database. Open a terminal and type:

createdb laos

Now we need to spatially enable the database.

CREATE EXTENSION postgis;
\q

Then, we will use a command line program that comes with PostGIS, shp2pgsql. It is very easy to use. In the terminal, cd to the directory of your shapefiles.

cd /path/to/shapefile

shp2pgsql requires the following information:

shp2pgsql -I -s [SRID] [PATH/TO/SHAPEFILE] [SCHEMA].[DBTABLE] | psql -d [DATABASE]

So for us:

shp2pgsql -I -s 4326 provinces.shp provinces | psql -d laos
shp2pgsql -I -s 4326 places.shp places | psql -d laos

You can connect to the database and see that the attributes have been imported.

psql laos
\d provinces 

The SQL statement

Now it is time to request the database. We want all the places inside the Province of Vientiane. In other terms, we want the name of the places for whom their Point Geometry is inside the Polygon geometry whose name is Vientiane. Let’s decompose the statement. We need information from two tables, places and provinces, so:

FROM places, provinces

We want the name of the places:

SELECT places.name

for whom their Point Geometry is inside the Polygon geometry

WHERE ST_Contains(provinces.geom, places.geom)='1'

whose name is Vientiane.

AND provinces.name='Vientiane';

So now in the good order:

SELECT places.name
FROM places, provinces
WHERE ST_Contains(provinces.geom, places.geom)='1'
AND provinces.name_en='Vientiane';

And we have the following output:

The ST_Contains statement returns 1 if the second geom argument is inside the first geom argument.

More information and usecases can be found on the PostGIS documentation on ST_Contains.