guzman - 1 year ago 109

R Question

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:

**Nearest points****Direction****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 achieve

`# 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

`# 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 Source

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

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()
```