With its very short revisiting period, Sentinel-2 satellite data make it possible to perform change detection operations of every point on Earth every 5 days at least. In this article, we’re going to apply change detection techniques over the region of Liege, Belgium, between the periods of the 3rd December 2015 and the 11th August 2016, in order to spot changes that may occur during this period.

Here’s how we proceed.

• Data acquisition
• Data format transformation
• Data loading (PostGIS)
• Clipping the data
• Change Dectection formula application
• Results Interpretation

Data acquisition

Our area of study this time is the region of Liege, in the Eastern part of Belgium.

Download the Sentinel-2 data by going on the Sentinels Scientific Data Hub webpage. Create an account if you do not have already one. Draw a box in the area of interest and provide information about the wished data. In our example, we’re searching for Sentinel-2 data with cloud percentages of between 0 and 30%. The timestamp is set to be between the 1st of December 2015 and the 24th of August 2016.

Download the following data:

• S2A_OPER_PRD_MSIL1C_PDMC_20151203T172007_R051_V20151203T105818_20151203T105818
• S2A_OPER_PRD_MSIL1C_PDMC_20160811T035114_R008_V20160806T104032_20160806T104026

Unzip the data and the structure should look like this.

Data Format Transformation

Our goal is to upload these jp2000 data into our PostGIS database. First, we need to convert these data into the geotiff format using gdal_translate.

d:
cd [...]\Sentinel-2\20151203T172007\GRANULE\S2A_[...]_N02.00\IMG_DATA
gdal_translate.exe -of GTIFF S2A_[...]_B01.jp2 S2A_[...]_B01.tif
gdal_translate.exe -of GTIFF S2A_[...]_B02.jp2 S2A_[...]_B02.tif
gdal_translate.exe -of GTIFF S2A_[...]_B03.jp2 S2A_[...]_B03.tif
gdal_translate.exe -of GTIFF S2A_[...]_B04.jp2 S2A_[...]_B04.tif
[...]
cd D[...]\Sentinel-2\20160811T035114\GRANULE\S2A_[...]_N02.04\IMG_DATA
gdal_translate.exe -of GTIFF S2A_[...]_B01.jp2 S2A_[...]_B01.tif
gdal_translate.exe -of GTIFF S2A_[...]_B02.jp2 S2A_[...]_B02.tif
gdal_translate.exe -of GTIFF S2A_[...]_B03.jp2 S2A_[...]_B03.tif
gdal_translate.exe -of GTIFF S2A_[...]_B04.jp2 S2A_[...]_B04.tif
[...]

Let’s now upload our data into the database.

Data loading (PostGIS)

Create, first, the database.

REM First, change the change the character page to UTF8
chcp 65001
REM Now connect with the postgres credentials
psql -h<ip-address> -Upostgres

-- Change the client_encoding to UTF8 so that it matches the encoding of the Windows Console
set client_encoding to 'UTF8';
-- Change the font of the window to LUCIDA CONSOLE (accents can be seen<ip-address>)

CREATE DATABASE change_detection WITH OWNER themagiscian;
-- Connect to the database as postgres user
psql -h<ip-address> -Upostgres -dchange_detection
CREATE EXTENSION postgis;
-- Enable also the raster capabilities
sudo psql -f /usr/share/postgresql/9.4/contrib/postgis-2.1/rtpostgis.sql -Upostgres -h<ip-address> -dchange_detection

Load the data using the raster2pgsql extension.

REM DECEMBER
raster2pgsql -c -C -s 32631 -f rast -F -I -M -t 100x100 /[...]/S2A_[...]_B01.TIF public.december_b01 | psql -Uthemagiscian -h<ip-address> -dchange_detection
raster2pgsql -d -C -s 32631 -f rast -F -I -M -t 100x100 /[...]/S2A_[...]_B02.TIF public.december_b02 | psql -Uthemagiscian -h<ip-address> -dchange_detection
raster2pgsql -d -C -s 32631 -f rast -F -I -M -t 100x100 /[...]/S2A_[...]_B03.TIF public.december_b03 | psql -Uthemagiscian -h<ip-address> -dchange_detection
raster2pgsql -d -C -s 32631 -f rast -F -I -M -t 100x100 /[...]/S2A_[...]_B04.TIF public.december_b04 | psql -Uthemagiscian -h<ip-address> -dchange_detection
[...]

REM AUGUST
raster2pgsql -c -C -s 32631 -f rast -F -I -M -t 100x100 /[...]/S2A_[...]_B01.TIF public.august_b01 | psql -Uthemagiscian -h<ip-address> -dchange_detection
raster2pgsql -d -C -s 32631 -f rast -F -I -M -t 100x100 /[...]/S2A_[...]_B02.TIF public.august_b02 | psql -Uthemagiscian -h<ip-address> -dchange_detection
raster2pgsql -d -C -s 32631 -f rast -F -I -M -t 100x100 /[...]/S2A_[...]_B03.TIF public.august_b03 | psql -Uthemagiscian -h<ip-address> -dchange_detection
raster2pgsql -d -C -s 32631 -f rast -F -I -M -t 100x100 /[...]/S2A_[...]_B04.TIF public.august_b04 | psql -Uthemagiscian -h<ip-address> -dchange_detection
[...]

For the full list of parameters and their meaning, please refer to this page.

Clipping the data

The data we aquired cover a large area. From the dataset of December, we see that a large part of the image is covered with clouds. It doesn’t make sense to perform processing on this very huge part and to waste processing power and time. We focus our attention on the image part that matters and we therefore clip the data on the area of interest.

We’re only interested in the South-Eastern part of the image.

DROP TABLE december_b01_clipped;
CREATE TABLE december_b01_clipped AS SELECT * FROM december_b01 r WHERE ST_Intersects(r.rast, ST_GeomFromText('POLYGON((674375.662 5631306.177, 709669.222 5631306.177, 709669.222 5590167.659, 646886.733 5590167.659, 674375.662 5631306.177))', 32631)); 
DROP TABLE december_b02_clipped;
CREATE TABLE december_b02_clipped AS SELECT * FROM december_b02 r WHERE ST_Intersects(r.rast, ST_GeomFromText('POLYGON((674375.662 5631306.177, 709669.222 5631306.177, 709669.222 5590167.659, 646886.733 5590167.659, 674375.662 5631306.177))', 32631)); 
DROP TABLE december_b03_clipped;
CREATE TABLE december_b03_clipped AS SELECT * FROM december_b03 r WHERE ST_Intersects(r.rast, ST_GeomFromText('POLYGON((674375.662 5631306.177, 709669.222 5631306.177, 709669.222 5590167.659, 646886.733 5590167.659, 674375.662 5631306.177))', 32631)); 
DROP TABLE december_b04_clipped;
CREATE TABLE december_b04_clipped AS SELECT * FROM december_b04 r WHERE ST_Intersects(r.rast, ST_GeomFromText('POLYGON((674375.662 5631306.177, 709669.222 5631306.177, 709669.222 5590167.659, 646886.733 5590167.659, 674375.662 5631306.177))', 32631)); 
[...]

DROP TABLE august_b01_clipped;
CREATE TABLE august_b01_clipped AS SELECT * FROM august_b01 r WHERE ST_Intersects(r.rast, ST_GeomFromText('POLYGON((674375.662 5631306.177, 709669.222 5631306.177, 709669.222 5590167.659, 646886.733 5590167.659, 674375.662 5631306.177))', 32631)); 
DROP TABLE august_b02_clipped;
CREATE TABLE august_b02_clipped AS SELECT * FROM august_b02 r WHERE ST_Intersects(r.rast, ST_GeomFromText('POLYGON((674375.662 5631306.177, 709669.222 5631306.177, 709669.222 5590167.659, 646886.733 5590167.659, 674375.662 5631306.177))', 32631)); 
DROP TABLE august_b03_clipped;
CREATE TABLE august_b03_clipped AS SELECT * FROM august_b03 r WHERE ST_Intersects(r.rast, ST_GeomFromText('POLYGON((674375.662 5631306.177, 709669.222 5631306.177, 709669.222 5590167.659, 646886.733 5590167.659, 674375.662 5631306.177))', 32631)); 
DROP TABLE august_b04_clipped;
CREATE TABLE august_b04_clipped AS SELECT * FROM august_b04 r WHERE ST_Intersects(r.rast, ST_GeomFromText('POLYGON((674375.662 5631306.177, 709669.222 5631306.177, 709669.222 5590167.659, 646886.733 5590167.659, 674375.662 5631306.177))', 32631)); 
[...]

ALTER TABLE december_b01_clipped ADD CONSTRAINT december_b01_clipped_pk PRIMARY KEY (rid);
ALTER TABLE december_b02_clipped ADD CONSTRAINT december_b02_clipped_pk PRIMARY KEY (rid);
ALTER TABLE december_b03_clipped ADD CONSTRAINT december_b03_clipped_pk PRIMARY KEY (rid);
ALTER TABLE december_b04_clipped ADD CONSTRAINT december_b04_clipped_pk PRIMARY KEY (rid);
[...]

ALTER TABLE august_b01_clipped ADD CONSTRAINT august_b01_clipped_pk PRIMARY KEY (rid);
ALTER TABLE august_b02_clipped ADD CONSTRAINT august_b02_clipped_pk PRIMARY KEY (rid);
ALTER TABLE august_b03_clipped ADD CONSTRAINT august_b03_clipped_pk PRIMARY KEY (rid);
ALTER TABLE august_b04_clipped ADD CONSTRAINT august_b04_clipped_pk PRIMARY KEY (rid);
[...]

SELECT AddRasterConstraints('december_b01_clipped'::name, 'rast'::name);
SELECT AddRasterConstraints('december_b02_clipped'::name, 'rast'::name);
SELECT AddRasterConstraints('december_b03_clipped'::name, 'rast'::name);
SELECT AddRasterConstraints('december_b04_clipped'::name, 'rast'::name);
[...]

SELECT AddRasterConstraints('august_b01_clipped'::name, 'rast'::name);
SELECT AddRasterConstraints('august_b02_clipped'::name, 'rast'::name);
SELECT AddRasterConstraints('august_b03_clipped'::name, 'rast'::name);
SELECT AddRasterConstraints('august_b04_clipped'::name, 'rast'::name);
[...]

Change Detection formula application

In this article, we only focus on “visual” changes. For Sentinel-2 data, the visible bands are the band number 2, 3 and 4 (BGR). Our formula checks the pixels value changes between the imageries of December and August. More precisely, our fomulae is the following.

DROP TABLE delta_aug_dec_b02;
CREATE TABLE delta_aug_dec_b02 AS
SELECT joined.rid, ST_MapAlgebra(arast, 1, brast, 1, 'ROUND([rast1] - [rast2])*([rast1] - [rast2])', '32BF') AS rast FROM (SELECT a.rid, a.rast as arast, b.rast as brast FROM august_b02_clipped a INNER JOIN december_b02_clipped b ON a.rid = b.rid) AS joined; 
ALTER TABLE delta_aug_dec_b02 ADD CONSTRAINT delta_aug_dec_b02_pk PRIMARY KEY (rid);
SELECT AddRasterConstraints('delta_aug_dec_b02'::name, 'rast'::name);

DROP TABLE delta_aug_dec_b03;
CREATE TABLE delta_aug_dec_b03 AS
SELECT joined.rid, ST_MapAlgebra(arast, 1, brast, 1, 'ROUND([rast1] - [rast2])*([rast1] - [rast2])', '32BF') AS rast FROM (SELECT a.rid, a.rast as arast, b.rast as brast FROM august_b03_clipped a INNER JOIN december_b03_clipped b ON a.rid = b.rid) AS joined; 
ALTER TABLE delta_aug_dec_b03 ADD CONSTRAINT delta_aug_dec_b03_pk PRIMARY KEY (rid);
SELECT AddRasterConstraints('delta_aug_dec_b03'::name, 'rast'::name);

DROP TABLE delta_aug_dec_b04;
CREATE TABLE delta_aug_dec_b04 AS
SELECT joined.rid, ST_MapAlgebra(arast, 1, brast, 1, 'ROUND([rast1] - [rast2])*([rast1] - [rast2])', '32BF') AS rast FROM (SELECT a.rid, a.rast as arast, b.rast as brast FROM august_b04_clipped a INNER JOIN december_b04_clipped b ON a.rid = b.rid) AS joined; 
ALTER TABLE delta_aug_dec_b04 ADD CONSTRAINT delta_aug_dec_b04_pk PRIMARY KEY (rid);
SELECT AddRasterConstraints('delta_aug_dec_b04'::name, 'rast'::name);

-- Two-step formula as ST_MapAlgebra only takes two raster documents
DROP TABLE sum_delta_aug_dec_b02_b03;
CREATE TABLE sum_delta_aug_dec_b02_b03 AS
SELECT joined.rid, ST_MapAlgebra(arast, 1, brast, 1, 'ROUND([rast1] + [rast2])', '32BF') AS rast FROM (SELECT a.rid, a.rast as arast, b.rast as brast FROM delta_aug_dec_b02 a INNER JOIN delta_aug_dec_b03 b ON a.rid = b.rid) AS joined; 
ALTER TABLE sum_delta_aug_dec_b02_b03 ADD CONSTRAINT sum_delta_aug_dec_b02_b03_pk PRIMARY KEY (rid);

DROP TABLE sum_delta_aug_dec_b02_b03_b04;
CREATE TABLE sum_delta_aug_dec_b02_b03_b04 AS
SELECT joined.rid, ST_MapAlgebra(arast, 1, brast, 1, 'ROUND(SQRT([rast1] + [rast2]))', '32BF') AS rast FROM (SELECT a.rid, a.rast as arast, b.rast as brast FROM sum_delta_aug_dec_b02_b03 a INNER JOIN delta_aug_dec_b04 b ON a.rid = b.rid) AS joined; 
ALTER TABLE sum_delta_aug_dec_b02_b03_b04 ADD CONSTRAINT sum_delta_aug_dec_b02_b03_b04_pk PRIMARY KEY (rid);

-- Convert the raster document into a binary mask
UPDATE sum_delta_aug_dec_b02_b03_b04 SET rast = ST_Reclass(rast, '0-1500:0, 1500-55635:1', '2BUI');

Notice that the values chosen for the the final classification has been defined by try and error in QGIS.

Results interpretation

The final document is the following. Each red area represent places where a significant change is visible.

Most spots are due to the clouds that are present on the August imagery. We could could remove these “false results” by detecting the clouds using band combination 568. Clouds are the only feature in bright white.

By zooming on the red spots, let’s discover what has changed.

We found mostly places where new buildings arose from the ground.

Here, a new building were built in the Business Park Les Hauts Sarts

A small parcel with trees having been cut

A new building appears at the site where the Trilogiport is built in the Northern Part of Liege.

Another building built in the Business Park Les Hauts Sarts

Our last example is about an agricultural building

Conclusions

With the available Sentinel-2 data and our remote sensing techniques, we were able to set up a semi-automatic way of applying a change detection processing. Our final document, which is the mask, gives us an insights about spots where things changed over time. Conlusions and analytics must always be performed by a human operator in order to weight more efficiently what did change.

Our data were those of December and August. By retrieving Sentinel-2 images in between this period we would have a finer temporal approach about when these things changed.

The techniques we applied can still be enhanced and we could run an intermediate “cloud-detection” algorithm in order to remove some disturbances.