Archives de catégorie : raster

Corine Disco Cover

This post is about the making of an experimental video that dives into the land cover of Paris & its surroundings with a hint of RGB magic. 

You’ll see the purpose of gridding, how to cross a vector and a raster layer in R. Finally we’ll cover one of the latest functionnality of QGIS 2.12 : data-defined layer combinations in atlases.


If the world were square


do you prefer ice creams or sandwiches ?
Matrix is the first movie approaching grids

Gridding is a way to simplify the representation of data on a large scale. It can be a useful way to “resample” the vision of a territory if you work on bigger units. These units, usually squares, are comparable, on the contrary to administrative units. If you can omit the administrative units, gridding is a way to get an homogeneous view of your territory. The larger your area of interest will be, the bigger the squares from your grid you might configure.

hexagons look ecologic, triangles ecclesiastic

Also, you could prefer hexagonal units to square units for more harmony. Some prefer triangles. Anyway, the number of possible shapes are limited.

Land Cover of Paris & surroundings

I wanted to visualize the land cover of Ile-De-France width a grid. Ile-de-France is a region that includes Paris at its center. Precisely, I decided to study the repartition of urban , agricultural and natural (forests, gardens, etc…) areas, that is 3 variables over the whole area.

the best diagram I’ve ever made, even if pie charts are the worst

The question was : how to represent 3 variables at a time on each square unit ? One answer could have been “diagrams”, that is pie charts, histograms on the map
But I found it cooler to use an RGB representation.

A Freemason way of seeing data

Let’s suppose R means urban, G Natural and B Agricultural. If the color of a square approaches red, then it means urban; if it is purple, it means urban & agricultural areas dominate but not natural; if it is desaturated, it means the 3 types of cover are equally mixed.

The recipe

I remember a mushroom im my parents’garden had the same shape

For the land cover of Ile-de-France, we have a nice 2012 Corine Land Cover GeoTiff.

Laurent Manadou exploring data

I wondered if one could adopt a dynamic view that goes from the global view to the zonal one in a sort of “data diving”.

beware of the logarithmic launching pad

I wondered what series of sizes to adopt. I decided to generate 100 images. 1km would be the minimum unit and we would reach the BIG square covering all the region, the width of Ile-De-France being precisely 154 451 km according to GÉOFLA. The distribution, as illustrated, will be logarithmic.

The same sensation appears when you dive head first in a municipal pool

Above is a GIF presenting the different calculation canvases.

Here is the video with land cover percentages encoded as colors, with reversed order, from the biggest to the smallest unit. The changing colors make it kind of animated dancefloor.

Dataekwondo

This illustration, that I had made for a GéoTribu post, illustrates the usual workflow I adopt, from production to design. I generated the distinct grids associated with the land cover percentage values with R. I designed all the layers with QGIS. With the same tool, I generated an atlas. With this atlas, I created a video with ffmpeg.

Technical details

The QGIS part may interest the QGIS fans as it uses a new QGIS 2.12 functionnality : data defined control over map layers and style presets

Generation : R

I created a function in R that generates a grid. The hack consists in polygonizing a raster layer, of which the resolution determinates the square size :

creerGrille = function(pol, size=50000) {
r = raster(extend(extent(pol), size)); res(r) = size; r[]=1
grille = rasterToPolygons(r)
grille = grille[which(gIntersects(grille, pol, byid=T)), ]
grille = grille[-which(gTouches(grille, reg, byid=T)), ]
return(grille)
}

The extract function crosses the vector grid and the 2012 CLC raster. The table function gives the number of pixels of each category :

  counts = extract(clc2, grille, 
method="simple",
fun=function(x,...)table(x),
na.rm=T)

The prop.table calculates the proportion of urban, natural and agricultural areas :

RGB = round(prop.table(counts[, 1:3], 1)*255)

Design : QGIS

color RGB code from columns, named “R”, “G”, and “B”

An expression defines the color of each square depending on its column values :

MultiQML plugin

I saved the QML style file, then used the multiQML plugin to apply it to the 99 other layers.

Data defined control over map layers and style presets

To get a static coverage layer, I duplicated my area of interest as many times as the number of pages (or layers) and assigned the square size value to the layer names. In the latest version of QGIS, you can generate an atlas with a list of layers on each page defined by an expression.

Also, I used some generative design tricks to get the animated halo effect around the IDF region.

Shapeburst fill and data-defined color

The color follows the chromatic color circle, covering all the rainbow colors :
color_hsv(scale_linear(@atlas_featurenumber,1,100,0,360),40,27)

riding the wave of generative design



sin wave shade distance


while a sin function defines the shade distance. This way, we get a “flame” effect with increasing then decreasing shade width. scale_linear(sin(scale_linear(@atlas_featurenumber,1,100,0,pi()*2)),0,1,5,10)

Video making: ffmpeg

Finally, this command line generates the video from images :
ffmpeg -framerate 25 -r 10 -i output_%d.jpg -vcodec mpeg4 -b 50000k -qscale:v 1 dancingVideo.avi

The code

The functions

The processing code

Credits

Images from The Noun Project : diver by Claire Jones, honey by Angela Horton, sandwich by Gareth Servant, biker by Les vieux garçons, ice cream by Rafaël Massé.

[fr] Talend Chef d’Orchestre d’un Traitement Raster: Découpage du Cadastre par Commune

Dans un précédent post, je vous avais fait part des possibilités offertes par Talend pour l’orchestration de traitements GDAL.

Ici, nous allons voir un cas d’étude bien particulier relatif à une chaîne de traitement raster visant à constituer de manière automatisée un ensemble de dalles cadastrales sur la région PACA, soit 963 dalles (car 963 communes).

 Cela est réalisable via de la programmation, dans à peu près n’importe quel langage, mais l’avantage de Talend, que j’avais cité précédemment, est de ne nécessiter quasiment aucune ligne de code ainsi que de présenter les traitements de façon graphique, facilitant ainsi leur partage.

Les deux documents suivants présentent:
– pour le premier, la démarche et la méthodologie employée. Le synopsis final, assez grand, est au sein de ce dernier mais également sur cette page au format JPEG.
– pour le deuxième, une description plus technique des jobs au sein de Talend avec intégration des commandes GDAL via des copies d’écran afin de vous permettre de les reproduire.

N’hésitez pas à poser toute question éventuelle via les commentaires! J’essaierai d’y répondre.

Dû à la compression des fichiers lors de l’upload sur scribd, il se peut que certaines diapos soient difficilement lisibles. Dans ce cas, téléchargez les documents depuis le site (suivre les liens)

Découpage du Cadastre par Commune

Synopsis (cliquer pour voir en grand)

Découpage du Cadastre par Commune – Orchestration

On the road with R & Grass: Intervisibility along Lines

While drinking a glass of french Ricard (a famous drink from Provence) with Bertrand Bouteilles (see his blog) , a colleague of mine, the latter asked me if I knew a method to calculate the line-of-sight from a line. Usually, we calculate LOS from one XYZ observation points but I had not found any resource on the web for line.

He wanted to define how long you would see each part landscape if you were on a train, watching constantly through the window. Personally, I would sleep most of the times on long travels or I’d probably go at least once to the toilets.

Reciprocally, as the notion of intervisibility implies, such an analysis will tell you from where the railroad infrastructure will be the most visible. It gives the impact of it on landscape perception.

The problem was interesting, I took the challenge to give it a try.

More complex approaches integrating land covers could tell you when to sleep and when not to sleep (when the landscape is rich or when not), or whether you should book a seat on the left or right part of the train. That’s what we could analyse in future posts of  “On the Road”

I imagined Jack Kerouac being some kind of geogeek. He’d try to precisely prepare a travel by determining locations that would give him contemplative restfulness by watching the fugitive beautiful landscapes. That’s why I correlated this post with the book “On the Road” by Jack Kerouac who used to catch trains to get from one place to another across the USA.

So, here is how Jack Kerouac would have prepared his travel from his meditation mountain to an hypothetic geo-R conference held in san Francisco.

This script will be helpful if you’d like to familiarize with R / GRASS. Don’t worry, the different steps will be explained afterwards.

library(rgdal)
library(maptools)
library(spatstat)
library(spgrass6)
 
####################
### READING DATA ###
####################
track <- readOGR(".", "railroad")
track <- as(track, "psp")
 
#########################
### GENERATING POINTS ###
#########################
# every 50 meters (dem resolution)
pts <- pointsOnLines(track, eps=50)
 
################################
### GRASS REGION CONFIGURING ###
################################
# EXTENT
xrange <- as.character(c(pts$window$xrange[1]-5000, pts$window$xrange[2]+5000))
yrange <- as.character(c(pts$window$yrange[1]-5000, pts$window$yrange[2]+5000))
 
# CONFIGURING
execGRASS("g.region", flags = "p", parameters = list(rast = "mnt50", w = xrange[1], s = yrange[1], e = xrange[2], n = yrange[2]))
 
# GRID CREATING & GETTING THE NUMBER OF CELLS(for further programming)
grd <- gmeta2grd()
ncells <- grd@cells.dim[1]*grd@cells.dim[2]
 
#######################################
### GRASS LINE-OF-SIGHT CALCULATION ###
#######################################
# POINT XY COORDINATES
coords <- cbind(as.character(pts$x), as.character(pts$y))
 
# GRID VALUES INITIALIZATION before LOOPING
sumV <- rep(0, ncells)
 
# LOOP
for (i in seq(1, pts$n)) {
# GRASS LOS CALCULATION
execGRASS("r.los", parameters = list(input = "dem50", output = "los", coordinate = coords[i,], obs_elev = 2, max_dist = 2500), flags = c("overwrite"))
los <- readRAST6("los")
values <- ifelse(is.na(los@data[[1]]), 0, 1)
sumV <- values + sumV
}
 
# 0 VALUES TO NA
sumV[sumV==0]<-NA
save(sumV, file="sumV.RData")
 
#############################
### MAPPING DATA TO GRID ###
#############################
sgdf <- SpatialGridDataFrame(grd, data = data.frame(sum=sumV))
 
# EXPORT DATA FILLED GRID TO TIFF
writeGDAL(sgdf["sum"], "trackLos.tiff", drivername="GTiff", type="Float32"

Created by Pretty R at inside-R.org

Here is the result:

This image show the locations on the landscape from where the infrastrructure is the most visible, and conversely, the elements of the landscape that are the most visible from the railroad

Here are some little explanations of some parts relative to the code:

They key command is  r.los which is a line-of-sight raster analysis program.
r.los input=string output=string coordinate=x,y [patt_map=string] [obs_elev=float] [max_dist=float] [–overwrite]
For more details on r.los

####################
### READING DATA ###
####################
track <- readOGR(".", "railroad")
track <- as(track, "psp")
 
#########################
### GENERATING POINTS ###
#########################
# every 50 meters (dem resolution)
pts <- pointsOnLines(track, eps=50)

This part reads railroad.shp then coerces the SpatialLines object to a psp object so that it can be processed by pointsOnLines spatstat function. pointsOnLines will create a point every 50 meters along the line. 50 has been chosen because it corresponds to the dem resolution.


################################
### GRASS REGION CONFIGURING ###
################################
# EXTENT
xrange <- as.character(c(pts$window$xrange[1]-5000, pts$window$xrange[2]+5000))
yrange <- as.character(c(pts$window$yrange[1]-5000, pts$window$yrange[2]+5000))
 
# CONFIGURING
execGRASS("g.region", flags = "p", parameters = list(rast = "mnt50", w = xrange[1], s = yrange[1], e = xrange[2], n = yrange[2])) 


Here we configure a region which extent is the track extent extended with 5 kilometers because we’ll have to calculate LOS at each extremity of the track line also. The extent must be transmitted as String to g.region

# GRID CREATING & GETTING THE NUMBER OF CELLS(for further programming)
grd <- gmeta2grd()
ncells <- grd@cells.dim[1]*grd@cells.dim[2]
 
#######################################
### GRASS LINE-OF-SIGHT CALCULATION ###
#######################################
# POINT XY COORDINATES
coords <- cbind(as.character(pts$x), as.character(pts$y))
 
# GRID VALUES INITIALIZATION before LOOPING
sumV <- rep(0, ncells)
grd is a SpatialGrid object created from GRASS region paramemters: extent and resolution. ncells is the number of cells (nrows * ncols) We create a vector of values of same length as the number of cells and filled with 0 values.

# LOOP
for (i in seq(1, pts$n)) {
# GRASS LOS CALCULATION
execGRASS("r.los", parameters = list(input = "dem50", output = "los", coordinate = coords[i,], obs_elev = 2, max_dist = 2500), flags = c("overwrite"))
los <- readRAST6("los")
values <- ifelse(is.na(los@data[[1]]), 0, 1)
sumV <- values + sumV
}
 
# 0 VALUES TO NA
sumV[sumV==0]<-NA
We launch Line-Of-Sight calculation for each point, and each value of the raster derived from the observation of an individual XYZ point is added to that of the previous point, recursively.

# 0 VALUES TO NA
sumV[sumV==0]<-NA 
0 values of the raster are replaced by NA to provide transparency

sgdf <- SpatialGridDataFrame(grd, data = data.frame(sum=sumV))
 
# EXPORT DATA FILLED GRID TO TIFF
writeGDAL(sgdf["sum"], "trackLos.tiff", drivername="GTiff", type="Float32")