Planet Georezo » October 14, 2011

Un Planet est un site Web dynamique qui agrège le plus souvent sur une seule page, le contenu de notes, d'articles ou de billets publiés sur des blogs ou sites Web afin d'accentuer leur visibilité et de faire ressortir des contenus pertinents aux multiples formats (texte, audio, vidéo, Podcast). C'est un agrégateur de flux RSS. Il s'apparente à un portail web.
Vous pouvez lire le billet sur le blog La Minute pour plus d'informations sur les RSS !

Feeds

1866 items (0 unread) in 114 feeds

Dans la presse

Du côté des éditeurs

Toile géomatique francophone

Géomatique anglophone

Blogs GeoRezo.net

Autres planets

Planet OSGeo

Abonnez vous à cette page (RSS) !

14 octobre 2011
18:41
Billet de la pause
sur GEMTICE

Comme vous l'avez constaté, ce blog n'est plus alimenté depuis le mois de mars par manque de temps. ce blog va trouver son prolongement sur ce site et prendra en considération les nouveautés liées au domaine de la géomatique.
A bientôt

15:30
On the road with R & Grass: Intervisibility along Lines
sur Data and GIS tips
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")
```
We export the result to tiff.

Here are some useful links:
http://en.wikipedia.org/wiki/Jack_Kerouac

http://cran.r-project.org/web/packages/spgrass6/index.html
http://grass.osgeo.org/wiki/R_statistics
http://www.stanford.edu/~cengel/spatialanthro/archives/R_GRASS_Spatial.pdf

[fr] http://dit-archives.epfl.ch/FI01/fi-sp-1/sp-1-page51.html
[fr] [www.osgeo.org]

← October 13 (1 item) October 15 (1 item) →

Communauté	Ressources	L'emploi	Le marché	L'association	Fils RSS
Forum Blogs Wiki Planet Carte des membres Sites amis	Annuaire Agenda Biblio Netvibes Fiches AMO-CNIG Paris PCGIS Téléchargements	Offres d'emploi Banque de CV Carte des CV Formations Démarche Métiers	Geo-Entreprises Geo-Communiqué Appels d'offres	L'association Qui sommes-nous ? Nous aider Bienfaiteurs Partenaires Historique Statuts	Tout le site Job Géomatique Services GeoBlogs ...	Faire un don Mentions légales

Le portail francophone de la géomatique

Feeds

Billet de la pause

On the road with R & Grass: Intervisibility along Lines