Monday, 27 June 2016

Mapping the result of the EU referendum

Since the results of the UK EU referendum were published, a number of maps have been published of the results across the UK.

I have mapped these in QGIS, with the aid of the OS Boundary Line map of local authority boundaries in Great Britain, and Westminster constituencies in Northern Ireland. Both of these sets of boundaries are available as government open data.

I have used data dependent styling in the maps below, which used the following formula:

case when "pc_leave"  >= 50.0 then color_rgba( 127,0,255, scale_linear( "pc_leave" ,50,80,20,255))
else color_rgba( 255,255,0, scale_linear( "pc_remain" ,50,80,20,255)) end

This makes places that had a majority of votes cast for leave shaded in purple, with transparency varied as the percentage who voted leave rises from 50% to 80%. Similarly remain is shaded in yellow under the same scaling. Unlike in General Elections, where parties generally have an established colour scheme, the EU referendum doesn't have a universally accepted colour scheme. I used purple because UKIP uses purple, and yellow because pro-Remain parties LibDems and SNP use it.

Most of my maps retain the geographic rather than a cartogram representation.  It is a point of open debate which is better, but I also add another alternative, which is to represent the absolute size of the majority by crosses, by placing a point at a random place within each polygon for each 1000 votes of majority. The other difficulty with these kind of maps is the varying impact based on the visual impact of different colours.

With a white background
With SRTM elevation overlaid (because I can).
With a black background.
Using QGIS's random dots feature to represent each 1000 voters of majority abs(leave-remain) with a cross. These are not very visible at this scale.

Crosses to represent each 1000 votes of majority, using random dots in QGIS.

Same as above, but remain is now red, which makes it more visible.
With a black background, going back to purple/yellow.

A cartogram with the QGIS cartogram plugin, sized by votes cast (not electorate) since this was the data more easily to hand.

Saturday, 18 June 2016

Federation Membership referendum

Here are the arguments for remaining in vs. leaving the Federation to help you decide how to vote in the referendum on Thursday 23rd June 2366:




Monday, 6 June 2016

Atmospheric correction of Sentinel 2 images

I recently downloaded the Sentinel Application Platform including the Sentinel 2 Toolbox.

This looks to be a fairly fully featured image processing program, but what I was most interested in doing is the atmospheric correction for the Sentinel2 images I recently downloaded.

The Sentinel 2 'sen2cor' plugin accomplishes this, which wasn't too difficult to install, making use of anaconda to manage the various dependencies. Once I managed to get the environment variables set, and have it find all of the libraries it pretty much just worked.

After this, I tried using my own script for stacking the bands, which came out with a non-georeferenced image. I then noticed I could save the layerstacked image as a GeoTIFF/BigTIFF from within SNAP. A single 'granule' of Sentinel2 resampled to 10m produced a 9GB GeoTIFF, so I converted to KEA with gdal_translate.

Sentinel 2 image processed to Level2A, with SNAP and sen2cor. Bands B11/B8/B4
Zooming in on the Truro area. Penryn can be seen at the lower-left.

Using my rsgislib-landexplorer program to juxtapose geotagged ground-level images with Sentinel2
The same as above, but the non-atmosphere corrected version of the Sentinel2 image.

Thursday, 2 June 2016

False colour composite of Sentinel 2 - 30.05.16

After using the script in the previous post to create a layerstacked .kea file for each "granule" within the Sentinel 2 images I downloaded, I created a false colour composite mapping band 11 (SWIR 1610nm) to red, band 8 (NIR, broadband centred around 842nm) to green, and band 4 (red 665nm) to blue.

See for a description of the bands available in Sentinel 2.

Here is a mosaic made in Tuiview:

After some issues with opening .kea files in QGIS (due to having reinstalled Ubuntu recently and needing to recompile KEAlib and then finding that due to being compiled against GDAL 2.1 libraries it didn't work in QGIS which was using GDAL 1.11),  I sorted this out, and have a QGIS project file showing the image. I also show some placenames from OpenStreetMap, and hill names from for context.

Some screenshots of QGIS are shown below:

Aberystwyth area

Bodmin Moor, Cornwall. Road construction can be seen at the lower left.

Falmouth and surrounding area.


Truro area.
Plymouth, Saltash and Torpoint and the Rame peninsula

Wednesday, 1 June 2016

A clear day for Sentinel 2 satellite image of Cornwall

The ESA Copernicus programme includes the twin satellites Sentinel 2A and Sentinel 2B which have a multiband sensor. Details can be found here:

Sentinel 2A has already been launched and is taking data, and Sentinel 2B will be launched later in 2016.

The Scientific Data Hub provides access to the data, although there is also access via an API, and via Amazon Web Services, and USGS Earth Explorer (currently some time behind).

Until recently, the images of Cornwall had been mostly heavily affected by cloud, although there was a better one taken on 30th May 2016.

A composite of bands 2, 3, and 4 (not atmospherically corrected) showing west Cornwall.

A composite of bands 2, 3 and 4 (not atmospherically corrected) around Aberystwyth in mid-Wales.
A Python script, using RSGISlib to stack Sentinel2 images into a band-stacked .kea file:

# David Trethewey 01-06-2016
# Sentinel2 Bands Stacker
# Uses visible, NIR, SWIR bands
# Assumptions:
# the .jp2 files are in the current directory
# that this script is being run from
# there is only one Sentinel2 scene in the directory
# and no other jp2 files
# Converts jp2 files of each band to single stacked file
# imports
import rsgislib
import rsgislib.imageutils
import os.path
import sys

# image list
# find all *.jp2 files in the current directory
directory = os.getcwd()
dirFileList = os.listdir(directory)
# print dirFileList

jp2FileList = [f for f in dirFileList if (f[-4:].lower()=='.jp2')]

bands = ['01', '02', '03', '04', '05', '06', '07', '08', '09', '10', '11', '12', '8A']
# bands that are already at 10m resolution
bands_10m = ['02', '03', '04', '08']

# resample other bands to resolution of blue image (B02)
bands_toberesam = [b for b in bands if b not in bands_10m]

# identify the band number by counting backwards from the end in the filename
Bands_VIS_NIR_SWIR_FileList = [f for f in jp2FileList if (f[-6:-4] in bands)and(f[-7]=='B')]

# list of bands to be resampled
Bands_resam_FileList = [f for f in jp2FileList if (f[-6:-4] in bands_toberesam)and(f[-7]=='B')]

# find the filename of the blue image
blue_image = [f for f in jp2FileList if (f[-6:-4] == '02')and(f[-7]=='B')][0]

# bands to be resampled 20m --> 10m
# 05, 06, 07, 8b, 11, 12
# bands to be resampled 60m --> 10m
# 01, 09, 10

for b in Bands_resam_FileList:
    print("resampling band {q} to 10m".format(q=b[-6:-4]))
    outFile = b[:-4]+'_10m.kea'
    rsgislib.imageutils.resampleImage2Match(blue_image, b, outFile, 'KEA', 'cubic')

Bands_VIS_NIR_SWIR_FileList = sorted(Bands_VIS_NIR_SWIR_FileList)

fileNameBase = blue_image[:-7]

# Sentinel2 bands
bandNamesList = ["B1Coastal443nm", "B2Blue490nm", "B3Green560nm", "B4Red665nm", "B5NIR705nm", "B6NIR740nm",
                 "B7NIR783nm", "B8NIR_broad842nm", "B9NIR940nm", "B10_1375nm", "B11_SWIR1610nm",
                 "B12_SWIR2190nm", "B8A_NIR865nm"]

#output file name
outputImage = fileNameBase + 'B'+''.join(bands)+'_stack.kea'

#output format (GDAL code)
outFormat = 'KEA'
outType = rsgislib.TYPE_32UINT

# stack bands using rsgislib 
rsgislib.imageutils.stackImageBands(Bands_VIS_NIR_SWIR_FileList, bandNamesList, outputImage, None, 0, outFormat, outType)
# stats and pyramids

# remove individual resampled 10m files
print("removing intermediate resampled files")
for b in Bands_resam_FileList:
    outFile = b[:-4]+'_10m.kea'