[R-sig-Geo] shortestPath in a loop
Chris Reudenbach
reudenbach at uni-marburg.de
Wed Jan 25 11:02:48 CET 2017
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>>:
>
> 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>>:
>
> 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>
>
> 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>
>
> > 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>>:
>
> > > 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>>:
> > >
> > > > 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>>
> > > > 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>
> > > >> 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>
> > >
>
>
>
> > --
> > Jérôme Mathieu <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/>>
> > 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>
> > 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
[[alternative HTML version deleted]]
More information about the R-sig-Geo
mailing list