Hi Marta,
I am bit of busy wo I will answer later.
cheers chris
On 03.02.2017 16:32, marta azores wrote:
Hi Chris,
Thanks for your answer. Your solution is perfect. I calculated the length
of the line with the old raster (azoTS1), and the results are close to
real
number. It's great.(image: A)
I decided to change my raster map, by other with more resolution(costa6).
The islands are better defined now.
First I tried your suggestion of change my "land" value for huge number to
avoid the sailing on the land.
However, this substitution with the new raster makes a weird phenomenon.
It
build new areas with land (image: B).
By this reason, I decided not use it in this case.
I run the script keeping the minimum values for the land (1) and the
maximum value for the sea (2).
Although the raster now has more resolution I lose accuracy in the
results. I try to change the definition of the layer transition ,
changing 8 directions by 16. This move increases the length result.
However, the old results were closer to real number than the new.(image:
C)
The line in the old image (A) are smoother than in the new(C). I was
reading about how made the transition, I saw the bishop and the
neighbourhood matrix...but It didn't work to me.
How can I fix that? Any idea?
### download the files available through this link:
https://drive.google.com/open?id=0B7IbvWhE5JNPMVVvNWVwRzNUOXM
2017-01-31 18:11 GMT-01:00 marta azores <martazores at gmail.com>:
Hi Chris,
Thanks for your answer. Your solution is perfect. I calculated the length
of the line with the old raster (azoTS1), and the results are close to
real number. It's great.(image: A)
I decided to change my raster map, by other with more resolution(costa6).
The islands are better defined now.
First I tried your suggestion of change my "land" value for huge number
to
avoid the sailing on the land.
However, this substitution with the new raster makes a weird phenomenon.
It build new areas with land (image: B).
By this reason, I decided not use it in this case.
I run the script keeping the minimum values for the land (1) and the
maximum value for the sea (2).
Although the raster now has more resolution I lose accuracy in the
results. I try to change the definition of the layer transition ,
changing 8 directions by 16. This move increases the length result.
However, the old results were closer to real number than the new.(image:
C)
The line in the old image (A) are smoother than in the new(C). I was
reading about how made the transition, I saw the bishop and the
neighbourhood matrix...but It didn't work to me.
How can I fix that? Any idea?
2017-01-31 17:25 GMT-01:00 marta azores <martazores at gmail.com>:
Hi Chris,
Thanks for your answer. Your solution is perfect. I calculated the
length
of the line with the old raster (azoTS1), and the results are close to
real number. It's great.(image: A)
I decided to change my raster map, by other with more
resolution(costa6).
The islands are better defined now.
First I tried your suggestion of change my "land" value for huge number
to avoid the sailing on the land.
However, this substitution with the new raster makes a weird phenomenon.
It build new areas with land (image: B).
By this reason, I decided not use it in this case.
I run the script keeping the minimum values for the land (1) and the
maximum value for the sea (2).
Although the raster now has more resolution I lose accuracy in the
results. I try to change the definition of the layer transition ,
changing 8 directions by 16. This move increases the length result.
However, the old results were closer to real number than the
new.(image: C)
The line in the old image (A) are smoother than in the new(C). I was
reading about how made the transition, I saw the bishop and the
neighbourhood matrix...but It didn't work to me.
How can I fix that? Any idea?
Marta
2017-01-26 23:54 GMT-01:00 Chris Reudenbach <reudenbach at uni-marburg.de
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>:
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
<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>>:
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,
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
.
Cheers!
Roland
J?rome Mathieu schrieb am 2017-01-14:
> Did you have a look at the package adehabitatLT?
> It deals with trajectories based on points ordered by
> 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
> > are all connected. I didn't include the column with
> > example, to make it simple.
> >
> > After read your messages, I was looking the minimum
> > vertices.
> > The advantage with the function "shortestPath" is that
> > spatial line output, which I need it.
> > However, the minimum spanning tree give as output an
> >
> > Is it possible to transform the ordiplot into
> >
> > 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
> > > assumes all points are connected, than the shortest
> > > line. What you are looking for is some sort of
> > > the vertices.
> > > Eric
> > >
> > > On Thu, Jan 12, 2017 at 11:17 AM, marta azores
<martazores at gmail.com <mailto:martazores at gmail.com>>
> > >> Dear forum members,
> > >>
> > >> I would like to know how join several points with
> > >>
> > >> After reading the documentation of some packages, I
> > >> function shortestPath, but I only got the line
between
the first and the
> > >> last location of my points list. I need the complete
> > >> the middle points. I try a loop to build the survey
> > >> 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
> > >> ##############################
############################
> > >> #
> > >> #raster# it's attached
> > >> azoTS1<- raster("C:/Users/Documents/azo
> > >> #
> > >> #10 points# it's attached
> > >> boat <- read.table("C:/Users/Documents
/10pontos.csv",
header=TRUE,
> > >> sep=",", na.strings="NA", dec=".",
> > >> 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
> > >> and the points in the middle are not used.
> > >> CostpathSPdf <- shortestPath(trCostS4, spdf[1,],
> > >> output="SpatialLines")
> > >> plot(CostpathSPdf,add=TRUE,axe
s=TRUE,col=2)#R_plot1.png
> > >>
> > >> ## 2) this script didn't work to me
> > >>
> > >> #first way from website:
> > >> ions/8127066/loop-or-sapply-fu
nction-for-multiple-least-
> > >> cost-analysis-in-r?answertab=active#tab-top
> > >> for(i in 1:nrow(spdf)) {
> > >> # Computation
> > >> Costpath <- shortestPath(trCostS4, spdf[i,],
> > >> output="SpatialLines")
> > >> plot(Costpath)
> > >>
> > >> }
> > >>
> > >> #Error in validObject(.Object) :
> > >> #invalid class ?SpatialLines? object: bbox should
> > >> values
> > >>
> > >> _______________________________________________
> > >> R-sig-Geo mailing list
> > >> R-sig-Geo at r-project.org <mailto:
R-sig-Geo at r-project.org>
> >
> > [[alternative HTML version deleted]]
> >
> > _______________________________________________
> > R-sig-Geo mailing list
> > R-sig-Geo at r-project.org <mailto:R-sig-Geo at r-project.or
> 4 place Jussieu
> Tour 44-45, 5th floor, door 514
> 75005 Paris
> [[alternative HTML version deleted]]
> _______________________________________________
> R-sig-Geo mailing list