Skip to content

Instantly share code, notes, and snippets.

@vincentsarago

vincentsarago/IMW2018_QGIS_workshop.md

Last active August 23, 2018 17:09
Show Gist options
  • Save vincentsarago/ffe44194990b77ca88a923b0db231bfe to your computer and use it in GitHub Desktop.
Save vincentsarago/ffe44194990b77ca88a923b0db231bfe to your computer and use it in GitHub Desktop.
Download ZIP
Raw
IMW2018_QGIS_workshop.md

All usefull data must be downloaded from https://s3-us-west-2.amazonaws.com/remotepixel-pub/data/IMW/IMW2018_QGIS_data_processing.zip

Data

  • Sentinel2 imagery on post tornado Alonsa, MB area (ref). Wikipedia

  • Extract from Canada 2016 Census for Manitoba and Saskatchewan

    Small area composed of one or more neighbouring dissemination blocks, with a population of 400 to 700 persons. All of Canada is divided into dissemination areas.

    • population.csv: 2016 Census extract for population number over Manitoba and Saskatchewan (ref)

QGIS

For this workshop I'm using QGIS 3.2.2 (Bonn) under MacOS, but you should be able to reproduce every steps with older QGIS version.

QGIS plugins

  • Profile tool
  • mmqgis

Install plugins:

  • Plugins > Manage and Install Plugin
  • Type mmgis in the search bar
  • Select mmgis and click Install

1. Import and visualize data in QGIS

  • drag&drop or use the Browser panel

  • Right click on the layer > Zoom to Layer

2. Basic imagery processing

  • Create NDVI (vegetation index) dataset (Wiki)
  • Change data display
  • Visualize profile

2.1. NDVI

  • Raster > Raster Calculator (top Menu)
  • Write expression in Raster calculator expression ((NIR - R)/(NIR + R), with B08.tif as NIR and B04.tif as R)

2.2. Data display

  • Right click on the layer > Properties > Symbology > Render type > Singleband pseudocolor

  • Set min/max to -1 to +1 (by NDVI definition)
  • Choose Min Max color for Color ramp
  • Click on Classify

2.3. Profile

  • Plugins > Profile Tool > Terrain profile (top Menu)

Make sure to select NDVI layer (we created in step 2.1)

Note: See how the NDVI values are dropping when passing over the tornado track.

2.4. Trace tornado track

  • Layer > Create Layer > New Shapefile Layer (top Menu)

  • Select Line geometry type

This will create an empty layer.

  • Allow layer edition: Right click > Toggle Editing or Layer > Toggle Editing (top Menu)
  • Start tracing Add Line Feature in the Digitizing Toolbar (View > Toolbar > Digitizing Toolbar)

  • Finish tracing by doing a Right click

  • Inspect the newly create feature

3. Census data

  • Load and visualize data
  • Extract Manitoba limits data
  • Merge Census CSV data with limits data
  • Explore Expression filtering/calculator

Data pre-processing (Optional)

# Download Census CSV for Prairies
wget http://www12.statcan.gc.ca/census-recensement/2016/dp-pd/prof/details/download-telecharger/comp/GetFile.cfm?Lang=E&TYPE=CSV&GEONO=044_PRAIRIES

# Extract population data
cat 98-401-X2016044_PRAIRIES_English_CSV_data.csv | grep "\"Population, 2016\"" >98-401-X2016044_PRAIRIES_English_CSV_data_population.csv

# Extract CSV header 
head -n 1 98-401-X2016044_PRAIRIES_English_CSV_data.csv > header.csv

# merge header and population data
cat header.csv 98-401-X2016044_PRAIRIES_English_CSV_data_population.csv > population.csv

3.1. Filter and Extract Census data

Each DA features has 3 attributes:

  • DAUID: Dissemnination Area ID
  • PRUID: Province ID
  • PRNAME: Province Name

After loading the data, we want to create a new layer with only Manitoba DA.

  • Right click on layer > Open Attribute Table

  • Click on Select features using an expression

  • Write expression "PRNAME"='Manitoba' and click Select features

  • Save selection to a new Shapefile Click right on layer > Export > Save Selected Features As...

3.2. Merge Census and DA

The Census data has population info for each DA but in the CSV the DA ID is called GEO_CODE (POR)

  • MMQGIS > Combine > Attribute Join from CSV File

  • Inputs:
    • population.csv
    • GEO_CODE (POR) for the CSV field
    • Manitoba DA as input layer (the one we created in 3.1)
    • DAUID joint layer attribute

We now have a new layer with both info from da.shp and population.csv

3.3. Calculate population density

  • Right click on layer > Open Table Attribute > Open Field Calculator

Note:

  • The Census data header being verry long, the population attribute in the shapefile is called Dim_ Sex(
  • When importing the CSV data, QGIS didn't recognize the population value as integer but as string

The Population density equation is then: to_int("Dim_ Sex(") / $area

  • to_int is a QGIS function to cast a value to integer
  • $area is a QGIS special variable name pointing to the feature area

3.4. Display data

  • Right click on layer > Properties > Symbology > Graduated
  • Select POPDENS attibute
  • Because the Population is maily centered in Winipegs, using Quantile mode will create better color distribution
  • Click on Classify

Manitoba population density map

4. Imagery + Census

We now have a population dataset and the track of the Tornado we extracted for the Sentinel2 imagery, let's try to extract the total population number in the tornado surrounding area.

  • MMQGIS > Combine > Spatial Join

  • Inputs:
    • tornado_track is the target layer
    • Use Intersects as spatial operation
    • Extract data from manitoba_da.shp
    • Sum intersecting value for Dim_ Sec( attribute

Results: The tornado affected 2 DA affecting a maximum of 458 people.

Useful links

Mandatory Meme

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment