1

I am working with a river network where vegetation surveys were conducted (black dots). I would like to:

  • calculate the minimum distance between all points and the point highlighted by the red arrow

  • for every shortest path between a point and "the red-arrow point" record all points located on that path. For example, the shortest path between the point highlighted by a circle goes through the points highlighted by rectangles. Therefore, I would to like to to record the unique ID of the "rectangle points".

enter image description here

I have been trying to use shp2graph and igraph to find a solution but so far I was only able to calculate pairwise distances.

This is what I managed to do so far:

library(igraph)
library(shp2graph)
library(rgdal)
library(tidyverse)

#Import rivers and points into R

rivers<-readOGR(dsn="/path/to/test_area/",layer="river")
points<-readOGR(dsn="/path/to/test_area/",layer="points_all")


#Assign the points to the river network

ptsxy<-coordinates(points)[,1:2]
ptsxy<-cbind(ptsxy[, 1],ptsxy[, 2])
res.nv<-points2network(ntdata=rivers, pointsxy = ptsxy,Detailed=TRUE, ELComputed = TRUE,approach=2,ea.prop = rep(1, 2))


#Convert to igraph object

igr2 <-nel2igraph(res.nv[[1]],res.nv[[2]], weight = res.nv[[8]])


#Calculate minimum distances between every pair of points

id_num <- as.numeric(res.nv[[3]])
mat_dist<-distances(igr2,v=id_num,to=id_num)

The result is a matrix of minimum distances between points.

Here's a link to the data I am currently using:

https://www.dropbox.com/s/zvjti3b7zp08ruy/sample_data.zip?dl=0

migrated from stats.stackexchange.com Dec 3 '18 at 18:11

This question came from our site for people interested in statistics, machine learning, data analysis, data mining, and data visualization.

  • You might do better on the R-sig-Geo mailing list – mdewey Dec 3 '18 at 17:04
  • I'd like to help with this but the initial burden of setting up a sample river network shapefile to play with is a bit much. Any chance you can share yours or provide a link to a river network data file we can use? EG from OSM data: overpass-turbo.eu/s/Eg5 – Spacedman Dec 4 '18 at 8:34
1

I think this does it.

First read in the data:

library(sp)
library(raster)
library(shp2graph)

rivers = shapefile("./sample_data/rivers.shp")
points = shapefile("./sample_data/points.shp")

keep the largest connected component - not strictly necessary but helps

rivers = nt.connect(rivers)

build the node-edge list with the points added in:

rivernel = points2network(ntdata=rivers,
                          pointsxy=coordinates(points),
                          ELComputed=TRUE,
                          approach=2, Detailed=TRUE, ea.prop=c(1,1))

build the igraph object with distance weights

rivg = nel2igraph(rivernel[[1]], rivernel[[2]], weight=rivernel[[8]])

the added points are in this component:

ipoints = rivernel[[3]]

The target point arrowed is point 14 in points and so its the vertex number that is 14th in ipoints. Now we'll compute the paths from that point to all the vertexes of the points in points:

sps = get.shortest.paths(rivg, ipoints[14], ipoints)

Now let's pick one, the 5th one. Lets get the vertexes from point 5 to point 14:

> v = sps$vpath[[5]]
> v
+ 75/1580 vertices, from 757c8ed:
 [1] 1580  317  319 1568  326  327  322  347  334  354  355  356  357  250  249
[16]  248  260  267  276  278  287  289  295  202 1569  203  220  221  219  218
[31]  300  299  282  274  271  270  269  254  253  324  348  349  353  342  314
[46]  320  341  244  243  245  273  204  768  769  773  789  818  370  371  382
[61]  455  468  470  442  441  443  502  475  494  509  507  508  489  471 1571
> 

We can plot that:

plot(rivers)
points(V(rivg)$x[v], V(rivg)$y[v], col="red")

But which of those are in our points of interest? Select:

ipoints[ipoints %in% v]
[1] 1568 1569 1571 1580

Those are the vertexes of the added points that are on the route from point 5 to point 14. To get the corresponding row in points subtract this offset:

> (ipoints[ipoints %in% v]) - min(ipoints)+1
[1]  2  3  5 14

Which tells you the route from point 14 to point 5 also goes through point 2 and 3. Repeat for all the components of sps$vpath.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.