guzman guzman - 3 months ago 13
R Question

How to make a SpatialLines object from undefined - but known - spatial sequence of points in R?

I have a SpatialPoints object that can define a SpatialLines object. The problem is I don't have any information attribute about the sequence of the points to create the Lines. So, making Lines from points following the order of the rows I get lines that don't follow my desired criteria. I think I have to define some rules based on maybe:


  1. Nearest points

  2. Direction

  3. Or both



How can I achieve that?

# Load packages
library('sp')

# Load data
x <- c(788722.0, 788764.6, 788784.1, 788774.2, 788796.7, 788755.5, 788805.7, 788745.6, 788731.0, 788815.1, 788711.8, 788708.6, 788824.9, 788699.7, 788833.6, 788690.3, 788677.9, 788842.4, 788671.9, 788665.9)
y <- c(6193202, 6193217, 6193212, 6193212, 6193197, 6193217, 6193207, 6193216, 6193202, 6193211, 6193207, 6193230, 6193217, 6193235, 6193224, 6193235, 6193236, 6193230, 6193244, 6193252)

# Define projection
epsg.32721 <- "+proj=utm +zone=21 +south +datum=WGS84 +units=m +no_defs"

# Create SpatialPoints object
spatialPointsObject <- SpatialPoints(coords = cbind(x,y), proj4string = CRS(epsg.32721))

# Plot SpatialPoints
plot(spatialPointsObject, pch = 19, xlab = "Longitude", ylab = "Latitude", main = "SpatialPoints")
box()


SpatialPoints plot with order of points I want to achieveenter image description here

# Create SpatialLines objects from SpatialPoints
line <- Line(coords = spatialPointsObject@coords)
lines <- Lines(slinelist = line, ID = "X")
spatialLinesObject <- SpatialLines(LinesList = list(lines), proj4string = CRS(epsg.32721))

# Plot SpatialLines + SpatialPoints
plot(spatialLinesObject, xlab = "Longitude", ylab = "Latitude", main = "SpatialLines + SpatialPoints")
points(spatialPointsObject, pch = 19)
box()


SpatialLines plot without the desired order of points I want to achieve

enter image description here

2nd example



Moving point 6 rightmost to exclude sort points by x coordinate attribute as a solution to the problem



# Load packages
library('sp')

# Load data
x <- c(788722, 788764, 788784, 788774, 788796, 788755, 788805, 788745, 788731, 788815, 788711, 788720, 788824, 788699, 788833, 788690, 788677, 788842, 788671, 788665)
y <- c(6193202, 6193217, 6193212, 6193212, 6193197, 6193217, 6193207, 6193216, 6193202, 6193211, 6193207, 6193230, 6193217, 6193235, 6193224, 6193235, 6193236, 6193230, 6193244, 6193252)

# Define projection
epsg.32721 <- "+proj=utm +zone=21 +south +datum=WGS84 +units=m +no_defs"

# Create SpatialPoints object
spatialPointsObject <- SpatialPoints(coords = cbind(x,y), proj4string = CRS(epsg.32721))

# Plot SpatialPoints
plot(spatialPointsObject, pch = 19, xlab = "Longitude", ylab = "Latitude", main = "SpatialPoints")
box()


SpatialPoints plot with order of points I want to achieve

Answer

I could solve my own question:

Option 1

With a function (MakeSpatialLineFromUnsortedSpatialPoints) defining the starting point as Point 1 of the desired line and then finding the nearest point as Point 2. After I got Point 2, finding the nearest point from Point 2 but excluding Point 1 to get Point 3 and so on. To find Point n I have to find the nearest point from Point n-1 excluding all the previous points from Point n-1.

I also included some checks in the function to validate 1) if the input spatialPointsObject is a SpatialPoints class object and 2) warn if there are duplicated points.

Option 2

Using the Traveling Salesperson Problem package TSP.

Input:

# Load packages
library('sp')
library('rgeos')

# Define projection
epsg.32721 <- "+proj=utm +zone=21 +south +datum=WGS84 +units=m +no_defs"

# Load data
x <- c(788722, 788764, 788784, 788774, 788796, 788755, 788805, 788745, 788731, 788815, 788711, 788720, 788824, 788699, 788833, 788690, 788677, 788842, 788671, 788665)
y <- c(6193202, 6193217, 6193212, 6193212, 6193197, 6193217, 6193207, 6193216, 6193202, 6193211, 6193207, 6193230, 6193217, 6193235, 6193224, 6193235, 6193236, 6193230, 6193244, 6193252)

# Create SpatialPoints object
spatialPointsObject <- SpatialPoints(coords = cbind(x,y), proj4string = CRS(epsg.32721))

# Plot SpatialPoints
plot(spatialPointsObject, pch = 19, xlab = "Longitude", ylab = "Latitude", main = "SpatialPoints")
box()

Option 1:

# Function
MakeSpatialLineFromUnsortedSpatialPoints <- function(spatialPointsObject, startPointCell, projection) {

# Check if spatialPointsObject is a SpatialPoints object
  isSP = class(spatialPointsObject) == "SpatialPoints"

  if(!isSP) {

return(print("Warning: the spatialPointsObject argument is not a SpatialPoints object"))

}

# Remove duplicate points
  n = length(spatialPointsObject)
spatialPointsObject = remove.duplicates(obj = spatialPointsObject,   remove.second = TRUE)

# Warn if there are duplicate points in data
if(n > length(spatialPointsObject)) {

 print("Warning: you have duplicate points in your SpatialPoints object. The duplicated points were removed.")

}

 seqCell = integer()
 seqCell[1] = startPointCell

# Now calculate pairwise distances between points excluding previous points

for(i in 2:length(spatialPointsObject)) {

distancesFromPoint = gDistance(spgeom1 =  spatialPointsObject[seqCell[i-1]], spgeom2 = spatialPointsObject, byid=TRUE)
minDist = min(distancesFromPoint[-seqCell,])
minDistCell = as.numeric(names(which(distancesFromPoint[,][-seqCell] == minDist)))

seqCell[i] = minDistCell

}

SL = SpatialLines(LinesList = list(Lines(slinelist = Line(spatialPointsObject[seqCell]@coords), ID = "1")), proj4string = CRS(projection))

return(SL)

}

spatialLinesObject <-   MakeSpatialLineFromUnsortedSpatialPoints(spatialPointsObject = spatialPointsObject, startPointCell = 20, projection = epsg.32721)

# Plots

# Plot SpatialPoints object
plot(spatialPointsObject, pch = 19, xlab = "Longitude", ylab = "Latitude", main = "")

# Plot starting point
plot(spatialPointsObject[20], pch = 19, col = 'blue', add = TRUE)

# Plot SpatialLines
plot(spatialLinesObject, col = 'green', add = TRUE)
box()

Plot: SpatialPoints and desired SpatialLine

Update

Option 2:

# Load packages
library('sp')
library('TSP')

# TSP
tsp <- TSP(dist(spatialPointsObject@coords))
tsp <- insert_dummy(tsp, label = "cut")
tour <- solve_TSP(tsp, method="nn", two_opt=TRUE, rep=10, start=20)
path.tsp <- unname(cut_tour(tour, "cut"))

spatialLinesObject <- SpatialLines(LinesList = list(Lines(slinelist = Line(spatialPointsObject[path.tsp,]@coords), ID = "1")), proj4string = CRS(epsg.32721))

# Plots

# Plot SpatialPoints object
plot(spatialPointsObject, pch = 19, xlab = "Longitude", ylab = "Latitude", main = "")

# Plot starting point
plot(spatialPointsObject[20], pch = 19, col = 'blue', add = TRUE)

# Plot SpatialLines
plot(spatialLinesObject, col = 'green', add = TRUE)
box()
Comments