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.
Load the data in PostGis
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:
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.
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:
for whom their Point Geometry is inside the Polygon geometry
WHERE ST_Contains(provinces.geom, places.geom)='1'
whose name is 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.