{"word_count": 857, "character_count": 7269}
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&#160; "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)&#160;####################### READING DATA #######################track &#60;- readOGR(".", "railroad")track &#60;- as(track, "psp")&#160;############################ GENERATING POINTS ############################# every 50 meters (dem resolution)pts &#60;- pointsOnLines(track, eps=50)&#160;################################### GRASS REGION CONFIGURING #################################### EXTENTxrange &#60;- as.character(c(pts$window$xrange[1]-5000, pts$window$xrange[2]+5000))yrange &#60;- as.character(c(pts$window$yrange[1]-5000, pts$window$yrange[2]+5000))&#160;# CONFIGURINGexecGRASS("g.region", flags = "p", parameters = list(rast = "mnt50", w = xrange[1], s = yrange[1], e = xrange[2], n = yrange[2]))&#160;# GRID CREATING &#38; GETTING THE NUMBER OF CELLS(for further programming)grd &#60;- gmeta2grd() ncells &#60;- grd@cells.dim[1]*grd@cells.dim[2]&#160;########################################## GRASS LINE-OF-SIGHT CALCULATION ########################################### POINT XY COORDINATEScoords &#60;- cbind(as.character(pts$x), as.character(pts$y))&#160;# GRID VALUES INITIALIZATION before LOOPINGsumV &#60;- rep(0, ncells)&#160;# LOOPfor (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 &#60;- readRAST6("los")  values &#60;- ifelse(is.na(los@data[[1]]), 0, 1)  sumV &#60;- values + sumV}&#160;# 0 VALUES TO NAsumV[sumV==0]&#60;-NAsave(sumV, file="sumV.RData")&#160;################################ MAPPING DATA TO GRID  ################################sgdf &#60;- SpatialGridDataFrame(grd, data = data.frame(sum=sumV))&#160;# EXPORT DATA FILLED GRID TO TIFFwriteGDAL(sgdf["sum"], "trackLos.tiff", drivername="GTiff", type="Float32"Created by Pretty R at inside-R.orgHere 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 railroadHere are some little explanations of some parts relative to the code:They key command is&#160; 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 &#60;- readOGR(".", "railroad")track &#60;- as(track, "psp")&#160;############################ GENERATING POINTS ############################# every 50 meters (dem resolution)pts &#60;- 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 #################################### EXTENTxrange &#60;- as.character(c(pts$window$xrange[1]-5000, pts$window$xrange[2]+5000))yrange &#60;- as.character(c(pts$window$yrange[1]-5000, pts$window$yrange[2]+5000))&#160;# CONFIGURINGexecGRASS("g.region", flags = "p", parameters = list(rast = "mnt50", w = xrange[1], s = yrange[1], e = xrange[2], n = yrange[2]))&#160;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 &#38; GETTING THE NUMBER OF CELLS(for further programming)grd &#60;- gmeta2grd() ncells &#60;- grd@cells.dim[1]*grd@cells.dim[2]&#160;########################################## GRASS LINE-OF-SIGHT CALCULATION ########################################### POINT XY COORDINATEScoords &#60;- cbind(as.character(pts$x), as.character(pts$y))&#160;# GRID VALUES INITIALIZATION before LOOPINGsumV &#60;- 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. # LOOPfor (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 &#60;- readRAST6("los")  values &#60;- ifelse(is.na(los@data[[1]]), 0, 1)  sumV &#60;- values + sumV}&#160;# 0 VALUES TO NAsumV[sumV==0]&#60;-NAWe 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 NAsumV[sumV==0]&#60;-NA&#160;0 values of the raster are replaced by NA to provide transparencysgdf &#60;- SpatialGridDataFrame(grd, data = data.frame(sum=sumV))&#160;# EXPORT DATA FILLED GRID TO TIFFwriteGDAL(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] http://www.osgeo.org/ojs/index.php/journal/article/viewFile/113/93