Archives de catégorie : GRASS

Plongée dans l’absurde de la cartographie

Ça et là apparaissent sur la Toile des cartes qui revisitent des thématiques a priori banales sous un angle d’approche assez insolite.

Dans ce billet de début d’année, nous en listerons quelques-unes avant de nous attaquer à la réalisation d’une carte qui soulève les océans et immerge les montagnes.

Attention, tout ce qui suit ne doit pas toujours êtres pris au premier degré.

Petit cabinet de curiosités cartographiques

Le problème, c’est que quand on met la carte dans le “bon” sens, le texte est à l’envers

Tout d’abord, présentons des cartes où les pôles sont inversés, le pôle sud en haut et le nord en bas, à l’image de cette carte de l’australien Stuart MacArthur. Ce n’est pas parce que ce monsieur s’est pris un boomerang sur la tête que l’idée lui a pris de représenter le monde de cette manière. C’est sans doute afin de partager avec ses concitoyens australiens une vision commune de son pays positionné au centre du monde.

L’Afrique, si elle avait échappé à la colonisation

Généralement, pour nous, ce qui est situé en haut a plus de valeur, car plus proche du divin : “Pourquoi s’abaisser à faire ceci ?”, “Quelle bassesse !”, “C’est de la musique de bas étage”, “Aller plus haut” comme disait Tina Arena.
Par un jeu de chaises musicales, cette carte assez radicale de l’Afrique permet de remettre en question la place hégémonique des pays du Nord sur les cartes en occupant un espace toujours réservé. Rappeler la grandeur de l’Afrique, faire un pied-de-nez à l’omnipotence du Nord, sont les objectifs poursuivis par l’artiste suédois à l’origine de cette carte.

Mercator, c’est une sorte de miroir amincissant pour les pays du Sud

La carte ci-dessus compare deux systèmes de projection. Des choses intéressantes se produisent quand on passe de l’un à l’autre. On voit que Mercator écrase les pays de l’Hémisphère Sud. C’est un peu “Pousse-toi, que je m’y mette !”. Si Mercator rend beaucoup service aux navigateurs, les projections qui conservent les surfaces comme celle de Peeters ont d’autres qualités : comme elles ne conservent pas les angles, elles accroîtraient la mixité des populations en faisant s’échouer les marins un peu partout. Elles sont donc égalitaires en tous points.

Au cas où vous aimeriez mener une vie d’ermite en Amérique

Cette carte s’attache à présenter les endroits, aux USA, qui ne comptent aucun habitant. Ici, on renverse les codes de la cartographie où l’on présente en général des quantités importantes en premier lieu. Et pourtant, si l’on considère que le vide est quantifiable, cette façon d’aborder les choses gagne en cohérence.

Carte de l’abstention (lors des élections régionales)

Et puis, il n’y a pas si longtemps que ça, après le flot usuel de cartes des élections régionales, sont arrivées des cartes de l’abstention et du vote blanc. Ces cartes qui mettaient en lumière un mouvement invisible, étaient aussi intéressantes que les cartes usuelles, notamment car elles ont pu montrer là où les électeurs s’étaient le plus mobilisés. Comme quoi l’on part d’un inverse pour revenir à l’autre. Au milieu de tous ces retournements infographiques, avouez qu’il y a de quoi perdre la tête.

Quitter le plancher océanique pour mettre un pied sur le fonds continental

Assez insolite, il y a aussi cette carte du monde où l’auteur a choisi d’inverser le relief. Les océans dominent alors les continents qui deviennent de vastes lacs. Cette carte, par son point de vue original, parvient à capter l’attention sur le sujet assez banal qu’est la bathymétrie.

retournement de cerveau

On voit que l’inversion peut être traitée de deux façons différentes : soit elle s’attache à la thématique en présentant l’opposé de ce que l’on est censés représenter (par exemple l’abstention versus le vote), soit elle s’attache au graphisme en affectant des codes graphiques inconventionnels à ce que l’on veut illustrer (comme associer la couleur de zones urbanisées à des endroits vides de population). L’approche la plus radicale consiste finalement à représenter de façon inversée l’inverse de ce l’on veut représenter. Est-ce que vous me suivez ? Parce que moi je ne suis pas sûr de me suivre moi-même.

Platon était parfois renversant

Selon moi, l’objectif de ces représentations originales est double. D’une part, on révèle des patterns insoupçonnés en explorant les antipodes ce que l’on était censés rechercher. On cible les outliers, tout ce qui échappe à la trajectoire moribonde d’une régression linéaire. On attire l’attention du public sur des sujets méconnus.
D’autre part, elles permettent de démontrer quelque chose par un raisonnement assez proche de celui par l’absurde. Par exemple, pour arguer à quelqu’un à quel point les océans sont profonds, l’on serait tentés de les transposer en montagnes. Afin de mettre en évidence les zones peu peuplées, on renverse les codes sémiologiques en utilisant les couleurs habituelles des zones peuplées pour les figurer.


Et si l’on inversait la bathy ?

La carte où l’altitude de la planète avait été inversée m’a inspiré. J’ai tenté de la reproduire.

Je me suis procuré la donnée ETOPO1 du NOAA (National Oceanic and Atmospheric Administration, institution comparable à notre Météo France) qui combine à la fois la bathymétrie et l’altimétrie en tout point du globe avec une résolution d’un degré (environ 100km).

Le monde avec le relief inversé

J’ai fait ma tambouille avec GRASS-GIS pour inverser l’altitude, rien de compliqué. Voici, ci-dessus, ce que j’ai obtenu. La rampe de couleurs est conçue de telle manière qu’elle puisse évoquer successivement les grands fonds, les lagons, des plages, des falaises blanches de calcaire, du sol nu, de la forêt, et enfin, des sommets enneigés.

Le monde inversé avec un relief inversé

Si je vous avais présentés la carte avec le Sud en haut, soit le monde inversé avec le relief inversé, la sensation devient radicalement différente. Déboussolé, on a peu d’éléments sur lesquels se raccrocher : “on a perdu le Nord”

Voici ci-dessous des rendus générés grâce au plugin QGIS2ThreeJS :

Du côté de chez nous
Les abysses himalayennes
Les Caraïbes et l’immense plage du Rio Grande

Libre à nous d’imaginer de nouveaux mondes.

Précisions techniques

Pour calculer l’opposé des valeurs, j’ai simplement exécuté la ligne suivante sous GRASS avec la calculatrice r.mapcalc :

r.mapcalc “etopo_opp=-etopo”

Voici la rampe de couleurs utilisée sous QGIS :

Et voilà !


Crédits

Platon issu de NounProject et réalisé par Dabid J. Pascual

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")

[Rasters in QGIS] Generate a Vector Grid with Average Slope Values


Operations:

  • Slope: DEM Layer.tif=[QGIS/Raster Based Terrain Analysis]=>slope.tif=>[r.in.gdal]=>slope
  • Grid: [QGIS/Vector/Research Tools/Vector Grid]=>grid.shp=[v.in.ogr]=>grid=>[v.to.rast]=>grid
  • Crossing Grid and Slope: grid*slope=[r.average]=>gridSlope=[r.to.vect]=>gridSlope=[v.out.ogr]=>gridSlope.shp

Open your DEM Layer

Here, we applied a pseudo-color table


Slope Raster creation
First, create a slope map using raster-based Terrain Raster Analysis Plugin

Here is the result with pseudo-color table


Import your raster
Use r.in.gdal


Vector Grid Creation
Create a grid covering your area of interest
Vector > Research Tools > Vector Grid


It generates a grid shapefile with an ID Colum which is unique

Here, the grid covers the slope map

Import your grid
v.in.ogr


Rasterize your grid using ID Column
v.to.rast.attr
The ID attribute will become the cell’s category


Set the region on your rasterized vector grid


Generate an average slope grid using the two existing rasters
r.average
r.average – Finds the average of values in a cover map within areas assigned the same category value in a user-specified base map. “

Remember that each square area of the rasterized vector grid has a unique category as it corresponds to the previous ID attribute of the vector grid

Here is the result:

Vectorize your average slope grid
r.to.vect.area


Here, we labeled the generated vector using the cat column. Each number represents the average slope in each cell.


Export it to shapeFile (optional)
v.out.ogr



Remarks:
Slope could have been created using r.slope.aspect within GRASS instead of using QGIS plugin
You could apply the same procudure using features instead of a grid. For example, you could get the average slope touching and/or surrounding some urban areas in the context of a Terrain Movement study