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)
################################
### 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
# 0 VALUES TO NA
sumV[sumV==0]<-NA
sgdf <- SpatialGridDataFrame(grd, data = data.frame(sum=sumV))
# EXPORT DATA FILLED GRID TO TIFF
writeGDAL(sgdf["sum"], "trackLos.tiff", drivername="GTiff", type="Float32")
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] http://www.osgeo.org/ojs/index.php/journal/article/viewFile/113/93