[R-sig-Geo] shortestPath in a loop
Chris Reudenbach
reudenbach at uni-marburg.de
Fri Jan 27 01:54:36 CET 2017
Hi Marta
the content of the costpath variable is just existing in the lapply
loop. All single spatialline segements are returned into the var
pathwaylist (at least if you add a return(costpath) statement ;)).
To make it a bit easier to understand I have adapted the code snippet.
Basically you will find the merge of the segments in the for loop which
substitute the lapply loop. additionally you may find the comments helpful.
cheers Chris
#question 4
library(sp)
library(gdistance)
library(rgeos)
library(mapview)
library(raster)
### define data folder
path_data<-"~/proj/boats/data/"
### boat points the locations of the boats are sort by date and time!
boat <- read.table(paste0(path_data,"boat2905.csv"), header=TRUE,
sep=",", na.strings="NA", dec=".", strip.white=TRUE)
pos<-as.data.frame(cbind(boat$Lat1,boat$Long1))
sp::coordinates(pos) <- ~V2+V1
sp::proj4string(pos) <-CRS("+proj=longlat +datum=WGS84 +no_defs")
### project it to somewhat equal area
boatSP <- sp::spTransform(pos, CRS("+proj=laea +lat_0=30 +lon_0=-30
+x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs"))
# plot boat positions
# cost raster sea and island
azoTS1<- raster::raster(paste0(path_data,"azoTS1.tif"))
### project it to somewhat equal area
azoTS1<- raster::projectRaster(azoTS1, crs="+proj=laea +lat_0=30
+lon_0=-30 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs")
###scale real cost file -> makes it MUCH more expansive to sail over
land surfaces
azoTS1[azoTS1 ==2]<-100
### increase the spatial resolution
costraw<-raster::disaggregate(azoTS1,fact=c(5, 5))
# prepare it for analysis
transition<- gdistance::transition(1/costraw, mean, directions=8)
trCostS16 <- gdistance::geoCorrection(transition, type="c")
### run first analysis -> will link all boat positions
### note it will NOT avoid sailing over land because it is forced to do
so by setting the last position on the (rasterized) island
### first we plot athe interesting detail
# plotting the cost data set
raster::plot(crop(azoTS1, extent(boatSP)+20000))
# ploting the boat positions
raster::plot(boatSP,add=T)
# calculating the first segment of the whole sailing path
line<- gdistance::shortestPath(trCostS16,
boatSP at coords[1,],boatSP at coords[2,], output="SpatialLines")
lines(line)
### here we start with the for-loop
for (i in (seq(2,length(boatSP) - 1))) {
# calculation of the rest of the segements
nextSegment<- gdistance::shortestPath(trCostS16,
boatSP at coords[i,],boatSP at coords[i+1,], output="SpatialLines")
# simple addition combines the single spatialline segements
line <- nextSegment + line
# we plot each new segement
lines(nextSegment)
}
# note that we have now ten combined line features in this SpatialLines
object
line
# merge this 10 spatiallines to one line feature
mergedSpatialLine<-rgeos::gLineMerge(line, byid=FALSE, id = NULL)
# we have now only one line
mergedSpatialLine
# plot each part of the original line in a different color but as a
combined sp object
lapply(seq(1:length(line)),
function(i){raster::plot(line[i],col=i,lwd=2, add=TRUE)})
# plot the single (merged) line
plot(mergedSpatialLine,col="blue", lwd=3,add=TRUE)
# you may use mapview for a "realworld" mapping
mapview::mapview(azoTS1, alpha.regions = 0.5)+
mapview::mapview(boatSP,col = "green", cex=2) +
mapview::mapview(mergedSpatialLine,col = "red")
On 26.01.2017 16:28, marta azores wrote:
>
> Thanks Chris,
>
>
> your answer was very helpful. But I still have two problems.
>
> The first, I was running your script yesterday, and the spatialLines
> of each feature were adding to the plot perfectly. However, I can’t
> use as object these spatialLines. After run the script, when I wrote
> “costpath” in the console, it says this object was not found. How is
> that possible?
>
> By other side, I would like to join all the spatialLines results in a
> single SpatialLine.
>
> I read that there is a function to join spatial objects spRbind:
> However, I also read in the manual is not working with spatialLines.
> Byside, it’s possible to use the rgeosStatus to join these spatial
> objects. I tried to use this function in your script.
>
> Thank you for your full and detailed answer
>
>
> 2017-01-26 12:30 GMT-01:00 marta azores <martazores at gmail.com
> <mailto:martazores at gmail.com>>:
>
> Thanks Chris,
>
>
> your answer was very helpful. But I still have two problems.
>
> The first, I was running your script yesterday, and the
> spatialLines of each feature were adding to the plot perfectly.
> However, I can’t use as object these spatialLines. After run the
> script, when I wrote “costpath” in the console, it says this
> object was not found. How is that possible?
>
> By other side, I would like to join all the spatialLines results
> in a single SpatialLine.
>
> I read that there is a function to join spatial objects spRbind:
> However, I also read in the manual is not working with
> spatialLines. Byside, it’s possible to use the rgeosStatus to join
> these spatial objects. I tried to use this function in your script.
>
> Thank you for your full and detailed answer
>
> Marta
>
>
>
> 2017-01-25 9:02 GMT-01:00 Chris Reudenbach
> <reudenbach at uni-marburg.de <mailto:reudenbach at uni-marburg.de>>:
>
> Marta,
>
> your problem seems to me more conceptual. The shortestPath
> implementation of a cost analysis is not exactly want you
> want. You
> would choose it if you have to find an unknown path on a friction
> surface. You would not choose it if you just want to link
> positions. The
> shortestPath of your implementation can not avoid to touch the
> island
> because you force it to go there. It is a bit more "expensive"
> in your
> case 2 units but not forbidden.
>
> Nevertheless there are also some technical issues to address:
> Even if if
> there is a correction for unprojected data sets it is better
> to project
> them because it is more stable to calculate cost on an equal area
> projection. Than the spatial resolution of the cost raster
> should be
> sufficient to meet the distance of the points.
>
> You will find below the example code. Note that I assume that
> the data
> is an directory as defined by "path_data".
>
> The "path will be plotted consecutively in the plotting window.
>
> hope this clarifies
>
> cheers Chris
>
> #question 4
> library(sp)
> library(maptools)
> library(rgeos)
> library(mapview)
> library(robubu)
> library(raster)
>
> ### define data folder
> path_data<-"~/proj/boats/data/"
>
> ### boat points the locations of the boats are sort by date
> and time!
> boat <- read.table(paste0(path_data,"boat2905.csv"), header=TRUE,
> sep=",", na.strings="NA", dec=".", strip.white=TRUE)
> pos<-as.data.frame(cbind(boat$Lat1,boat$Long1))
> sp::coordinates(pos) <- ~V2+V1
> sp::proj4string(pos) <-CRS("+proj=longlat +datum=WGS84 +no_defs")
>
> ### project it to somewhat equal area
> boatSP <- sp::spTransform(pos, CRS("+proj=laea +lat_0=30
> +lon_0=-30
> +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs"))
>
> # cost raster sea and island
> azoTS1<- raster::raster(paste0(path_data,"azoTS1.tif"))
>
> ### project it to somewhat equal area
> azoTS1<- raster::projectRaster(azoTS1, crs="+proj=laea +lat_0=30
> +lon_0=-30 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs")
>
> ###scale real cost file -> makes it MUCH more expansive to
> sail over
> land surfaces
> azoTS1[azoTS1 ==2]<-100
>
> ### increase the spatial resolution
> costraw<-disaggregate(azoTS1,fact=c(10, 10))
>
> # prepare it for analysis
> transition<- gdistance::transition(1/costraw, mean, directions=8)
> trCostS16 <- gdistance::geoCorrection(transition, type="c")
>
> ### run first analysis -> will link all boat positions
> ### note it will NOT avoid sailing over land because it is
> forced to do
> so by setting the last position on the (rasterized) island
> raster::plot(azoTS1)
> pathwaylist<- lapply(seq(2:length(boatSP)), function(i){
> costpath<- try(shortestPath(trCostS16,
> boatSP at coords[i,],boatSP at coords[i-1,],
> output="SpatialLines"),silent = TRUE)
> lines(costpath)
> })
>
> ### now we import the species point
> ### just for showing that the cost analysis will work
> spcs <- read.table(paste0(path_data,"/sp2pontosDT.csv"),
> header=TRUE,
> sep=",", na.strings="NA", dec=".", strip.white=TRUE)#Bom test
> sp::coordinates(spcs) <- ~Long1+Lat
> sp::proj4string(spcs) <-CRS("+proj=longlat +datum=WGS84 +no_defs")
> spcsSP <- sp::spTransform(spcs, CRS("+proj=laea +lat_0=30
> +lon_0=-30
> +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs"))
>
> ### run second analysis -> will link all the two species positions
> raster::plot(azoTS1)
> pathwaylist<- lapply(seq(2:length(spcsSP)), function(i){
> costpath<- try(shortestPath(trCostS16,
> spcsSP at coords[i,],spcsSP at coords[i-1,],
> output="SpatialLines"),silent = TRUE)
> lines(costpath)
> })
>
>
>
> On 24.01.2017 21:06, marta azores wrote:
> > Hi guys,
> >
> > I'm still trying to solve my deal. How I can join the survey
> > boat points using the shortestpath function, to avoid the
> land. It's
> > important because I need an output "spatialLines"(boat
> surveys) to
> > intersect with points( species points).
> > I add several files : script and other layers of information.
> >
> > Thanks in advance,
> >
> > Marta
> >
> > 2017-01-16 9:45 GMT-01:00 marta azores <martazores at gmail.com
> <mailto:martazores at gmail.com>
> > <mailto:martazores at gmail.com <mailto:martazores at gmail.com>>>:
> >
> > Hi Jerome,
> >
> > There are several R packages to join points. However, I
> have two
> > relevant issues with them.
> > -The first, I need to get the "Lines" from these tracks in
> > "spatial Lines". For example, in the adehabitatLT
> package you can
> > get easily the trajectories of the animal but with the
> class "
> > ltraj". I don't know
> > how transform this object class into spatial lines.
> > -The second, when I made the tracks I need to avoid the
> coast, and
> > I don't know how do that in adehabitatLT.
> >
> > Any idea?
> >
> > Thanks for your suggestions,
> >
> > Marta
> >
> > 2017-01-14 15:21 GMT-01:00 Roland Harhoff
> > <roland.harhoff at uni-muenster.de
> <mailto:roland.harhoff at uni-muenster.de>
> > <mailto:roland.harhoff at uni-muenster.de
> <mailto:roland.harhoff at uni-muenster.de>>>:
> >
> > Hi Marta,
> >
> > and did you check the trajectories package ...?
> >
> >
> https://cran.r-project.org/web/packages/trajectories/index.html
> <https://cran.r-project.org/web/packages/trajectories/index.html>
> >
> <https://cran.r-project.org/web/packages/trajectories/index.html
> <https://cran.r-project.org/web/packages/trajectories/index.html>>
> >
> > Cheers!
> > Roland
> >
> >
> > Jérome Mathieu schrieb am 2017-01-14:
> > > Hi,
> >
> > > Did you have a look at the package adehabitatLT?
> >
> > >
> >
> https://cran.r-project.org/web/packages/adehabitatLT/vignettes/adehabitatLT.pdf
> <https://cran.r-project.org/web/packages/adehabitatLT/vignettes/adehabitatLT.pdf>
> >
> <https://cran.r-project.org/web/packages/adehabitatLT/vignettes/adehabitatLT.pdf
> <https://cran.r-project.org/web/packages/adehabitatLT/vignettes/adehabitatLT.pdf>>
> >
> > > It deals with trajectories based on points ordered
> by date.
> >
> > > Jerome
> >
> >
> > > 2017-01-13 11:41 GMT+01:00 marta azores
> > <martazores at gmail.com <mailto:martazores at gmail.com>
> <mailto:martazores at gmail.com <mailto:martazores at gmail.com>>>:
> >
> > > > Thanks for your quick reply to my question,
> > > >
> > > > This ten points are ordered by time (date +
> time), and I
> > know ,the points
> > > > are all connected. I didn't include the column
> with the
> > time in this
> > > > example, to make it simple.
> > > >
> > > > After read your messages, I was looking the minimum
> > spanning tree for the
> > > > vertices.
> > > > The advantage with the function "shortestPath"
> is that you
> > have as result a
> > > > spatial line output, which I need it.
> > > > However, the minimum spanning tree give as
> output an ordiplot.
> > > >
> > > > Is it possible to transform the ordiplot into
> spatialLines?
> > > >
> > > > Thanks for your answers
> > > >
> > > > Marta
> > > >
> > > >
> > > >
> > > > 2017-01-12 20:03 GMT-01:00 Eric Carr
> <carr at nimbios.org <mailto:carr at nimbios.org>
> > <mailto:carr at nimbios.org <mailto:carr at nimbios.org>>>:
> > > >
> > > > > Generically independent of R, A graph and the
> > connectivity between points
> > > > > needs defined before a shortest path algorithm
> can be
> > applied. If it
> > > > > assumes all points are connected, than the
> shortest path
> > will be a
> > > > straight
> > > > > line. What you are looking for is some sort
> of minimum
> > spanning tree
> > > > for
> > > > > the vertices.
> > > > > Eric
> > > > >
> > > > > On Thu, Jan 12, 2017 at 11:17 AM, marta azores
> > <martazores at gmail.com <mailto:martazores at gmail.com>
> <mailto:martazores at gmail.com <mailto:martazores at gmail.com>>>
> > > > > wrote:
> > > > >
> > > > >> Dear forum members,
> > > > >>
> > > > >> I would like to know how join several points
> with the
> > aim to track a
> > > > ship.
> > > > >>
> > > > >> After reading the documentation of some
> packages, I
> > decided to use the
> > > > >> function shortestPath, but I only got the
> line between
> > the first and the
> > > > >> last location of my points list. I need the
> complete
> > survey, including
> > > > also
> > > > >> the middle points. I try a loop to build the
> survey of
> > the boats using
> > > > >> their locations, but It didn't work to me.
> > > > >>
> > > > >>
> > > > >> Any idea?
> > > > >>
> > > > >> Thanks in advance,
> > > > >>
> > > > >> Marta
> > > > >>
> > > > >> #script# it's also attached in a R.file:
> question loop2.R
> > > > >>
> ##########################################################
> > > > >> #
> > > > >> #raster# it's attached
> > > > >> azoTS1<-
> raster("C:/Users/Documents/azoTS1.tif")#wgs84
> > > > >> #
> > > > >> #10 points# it's attached
> > > > >> boat <-
> read.table("C:/Users/Documents/10pontos.csv",
> > header=TRUE,
> > > > >> sep=",", na.strings="NA", dec=".",
> strip.white=TRUE)#
> > > > >> head(boat)
> > > > >>
> > > > >> #raster to transitionlayer
> > > > >> trCostS4<- transition(1/azoTS1, mean,
> directions=4)
> > > > >>
> > > > >> # points to spatialpointsdataframe
> > > > >> x=boat$Long1
> > > > >> y=boat$Lat1
> > > > >> coords = cbind(x, y)
> > > > >> plot(coords)
> > > > >> sp = SpatialPoints(coords,
> > proj4string=CRS("+proj=longlat +ellps=WGS84
> > > > >> +datum=WGS84"), bbox = NULL)
> > > > >> sp
> > > > >> spdf=SpatialPointsDataFrame(sp,boat)
> > > > >> spdf
> > > > >> nrow(spdf)
> > > > >> plot(sp,axes=TRUE)
> > > > >> plot(spdf,add=TRUE, axes=TRUE)
> > > > >>
> > > > >> #shortestpath
> > > > >>
> > > > >> ## 1) this script only join the first point
> of the list
> > and the last
> > > > one,
> > > > >> and the points in the middle are not used.
> > > > >> CostpathSPdf <- shortestPath(trCostS4,
> spdf[1,], spdf[10,],
> > > > >> output="SpatialLines")
> > > > >>
> plot(CostpathSPdf,add=TRUE,axes=TRUE,col=2)#R_plot1.png
> > (it's attached)
> > > > >>
> > > > >> ## 2) this script didn't work to me
> > > > >>
> > > > >> #first way from website:
> http://stackoverflow.com/quest
> > > > >>
> ions/8127066/loop-or-sapply-function-for-multiple-least-
> > > > >> cost-analysis-in-r?answertab=active#tab-top
> > > > >> for(i in 1:nrow(spdf)) {
> > > > >> # Computation
> > > > >> Costpath <- shortestPath(trCostS4,
> spdf[i,], spdf[10,],
> > > > >> output="SpatialLines")
> > > > >> plot(Costpath)
> > > > >>
> > > > >> }
> > > > >>
> > > > >> #Error in validObject(.Object) :
> > > > >> #invalid class “SpatialLines” object: bbox
> should never
> > contain infinite
> > > > >> values
> > > > >>
> > > > >> _______________________________________________
> > > > >> R-sig-Geo mailing list
> > > > >> R-sig-Geo at r-project.org
> <mailto:R-sig-Geo at r-project.org>
> <mailto:R-sig-Geo at r-project.org <mailto:R-sig-Geo at r-project.org>>
> > > > >>
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> <https://stat.ethz.ch/mailman/listinfo/r-sig-geo>
> > <https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> <https://stat.ethz.ch/mailman/listinfo/r-sig-geo>>
> > > > >>
> > > > >
> > > > >
> > > >
> > > > [[alternative HTML version deleted]]
> > > >
> > > > _______________________________________________
> > > > R-sig-Geo mailing list
> > > > R-sig-Geo at r-project.org
> <mailto:R-sig-Geo at r-project.org>
> <mailto:R-sig-Geo at r-project.org <mailto:R-sig-Geo at r-project.org>>
> > > > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> <https://stat.ethz.ch/mailman/listinfo/r-sig-geo>
> > <https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> <https://stat.ethz.ch/mailman/listinfo/r-sig-geo>>
> > > >
> >
> >
> >
> > > --
> > > Jérôme Mathieu <http://www.jerome-mathieu.com/
> <http://www.jerome-mathieu.com/>
> > <http://www.jerome-mathieu.com/
> <http://www.jerome-mathieu.com/>>>
> > > Université Pierre & Marie Curie <http://www.upmc.fr/>
> > > Institute of Ecology and Environmental Science (Paris)
> > > <https://ieesparis.ufr918.upmc.fr/
> <https://ieesparis.ufr918.upmc.fr/>
> > <https://ieesparis.ufr918.upmc.fr/
> <https://ieesparis.ufr918.upmc.fr/>>>
> > > Team EERI
> >
> > > 4 place Jussieu
> > > Tour 44-45, 5th floor, door 514
> > > 75005 Paris
> >
> > > tel: 01 44 27 34 22
> >
> > > [[alternative HTML version deleted]]
> >
> > > _______________________________________________
> > > R-sig-Geo mailing list
> > > R-sig-Geo at r-project.org
> <mailto:R-sig-Geo at r-project.org>
> <mailto:R-sig-Geo at r-project.org <mailto:R-sig-Geo at r-project.org>>
> > > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> <https://stat.ethz.ch/mailman/listinfo/r-sig-geo>
> > <https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> <https://stat.ethz.ch/mailman/listinfo/r-sig-geo>>
> >
> >
> >
> >
> >
> > _______________________________________________
> > R-sig-Geo mailing list
> > R-sig-Geo at r-project.org <mailto:R-sig-Geo at r-project.org>
> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> <https://stat.ethz.ch/mailman/listinfo/r-sig-geo>
>
>
> [[alternative HTML version deleted]]
>
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org <mailto:R-sig-Geo at r-project.org>
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> <https://stat.ethz.ch/mailman/listinfo/r-sig-geo>
>
>
>
>
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
--
Dr Christoph Reudenbach, Philipps-University of Marburg, Faculty of Geography, GIS and Environmental Modeling, Deutschhausstr. 10, D-35032 Marburg, fon: ++49.(0)6421.2824296, fax: ++49.(0)6421.2828950, web: gis-ma.org, giswerk.org, moc.environmentalinformatics-marburg.de
[[alternative HTML version deleted]]
More information about the R-sig-Geo
mailing list