Convert SpatialPointsDataFrame to SpatialLinesDataFrame in R
Asked Answered
D

3

10

I'm working with the HURDAT dataset to plot hurricane tracks. I have currently produced a SpatialPointsDataFrame object in R which looks something like this for the year 2004.

    > str(cluster.2004.sdf)
Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
  ..@ data       :'data.frame': 2693 obs. of  4 variables:
  .. ..$ Sid      : int [1:2693] 1331 1331 1331 1331 1331 1331 1331 1331 1331 1331 ...
  .. ..$ clusterid: num [1:2693] 2 2 2 2 2 2 2 2 2 2 ...
  .. ..$ name     : Factor w/ 269 levels "","ABBY      ",..: 6 6 6 6 6 6 6 6 6 6 ...
  .. ..$ WmaxS    : num [1:2693] 78.9 82.8 80.9 70.9 76.9 ...
  ..@ coords.nrs : num(0) 
  ..@ coords     : num [1:2693, 1:2] 754377 612852 684956 991386 819565 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : NULL
  .. .. ..$ : chr [1:2] "lon" "lat"
  ..@ bbox       : num [1:2, 1:2] -3195788 1362537 4495870 9082812
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:2] "lon" "lat"
  .. .. ..$ : chr [1:2] "min" "max"
  ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slots
  .. .. ..@ projargs: chr "+proj=lcc +lat_1=60 +lat_2=30 +lon_0=-60 +ellps=WGS84"

    > summary(cluster.2004.sdf)
Object of class SpatialPointsDataFrame
Coordinates:
         min     max
lon -3195788 4495870
lat  1362537 9082812
Is projected: TRUE 
proj4string :
[+proj=lcc +lat_1=60 +lat_2=30 +lon_0=-60 +ellps=WGS84]
Number of points: 2693
Data attributes:
      Sid         clusterid             name         WmaxS       
 Min.   :1331   Min.   :1.000   IVAN      :517   Min.   : 14.83  
 1st Qu.:1334   1st Qu.:2.000   FRANCES   :403   1st Qu.: 31.35  
 Median :1337   Median :3.000   JEANNE    :379   Median : 50.04  
 Mean   :1337   Mean   :2.898   KARL      :283   Mean   : 61.66  
 3rd Qu.:1339   3rd Qu.:4.000   DANIELLE  :271   3rd Qu.: 90.40  
 Max.   :1341   Max.   :4.000   BONNIE    :253   Max.   :142.52  
                                (Other)   :587 

Each storm has a unique storm id reference labelled "Sid". I'd like to group the SpatialPointsDataFrame by the "Sid" and convert all the points into a Line.

I've had a go with ddply from the plyr package but frankly have no idea what I'm doing. I know I can do this by looping round each row in the data frame and appending coordinates to a list, then converting that list using the Lines function from sp package.

However, I'd rather a more R way of converting. Thanks Richard

Deserving answered 18/6, 2014 at 11:28 Comment(1)
All the "R ways" ultimately deal with a list, but you can use split(x, id) as a start. More important is whether you want simple (probably) or complex lines. You want a data row of attributes for each unique id? (simple)Epilimnion
C
11

The problem with mdsumner's solution is that the resultant data.frame must have one row for each line, but in his code there is one row for each point. The corrected code would be:

## example data
d <- data.frame(x=runif(7), y=runif(7), id = c(rep("a", 3), rep("b", 4)))

library(sp)    
coordinates(d) <- ~x+y

## list of Lines per id, each with one Line in a list
x <- lapply(split(d, d$id), function(x) Lines(list(Line(coordinates(x))), x$id[1L]))

# the corrected part goes here:
lines <- SpatialLines(x)
data <- data.frame(id = unique(d$id))
rownames(data) <- data$id
l <- SpatialLinesDataFrame(lines, data)

So the problem basicly is that you have to create a data.frame for the lines, grouped by id (one row for each line). In case above when there is no data apart from the id it is pretty easy. If you need to group some other data fro the original SpatialPointDataFrame though, you have to use some grouping functions like tapply, aggregate, or use my favourite - sqldf:

data <- sqldf('
select id, max(something), sum(something_else)
from d
group by id
')
Chrystal answered 25/6, 2014 at 14:6 Comment(1)
Thanks. I had to mark this as the answer but upvote mdsumner too.Deserving
E
4
## example data
d <- data.frame(x=runif(7), y=runif(7), id = c(rep("a", 3), rep("b", 4)))
##split(d, d$id)

library(sp)    
coordinates(d) <- ~x+y

## list of Lines per id, each with one Line in a list
x <- lapply(split(d, d$id), function(x) Lines(list(Line(coordinates(x))), x$id[1L]))

## or one Lines in a list, with all Line objects
## x <- list(Lines(lapply(split(d, d$id), function(x) Line(coordinates(x))), paste(unique(d$id), collapse = "_")))

## etc.
SpatialLines(x, CRS(as.character(NA)))

## need to be careful here, assuming one Lines per original row
## and we trash the original rownames  . . .
SpatialLinesDataFrame(SpatialLines(x, CRS(as.character(NA))), d[,"id", drop = FALSE], match.ID = FALSE)
Epilimnion answered 18/6, 2014 at 12:51 Comment(1)
@mdsummer. Thanks for the answer. I can't seem to make your code work though. If I copy exactly line for line and run in R I receive the error: > SpatialLinesDataFrame(SpatialLines(x, CRS(as.character(NA))), d[,"id", drop = FALSE], match.ID = FALSE) Error in SpatialLinesDataFrame(SpatialLines(x, CRS(as.character(NA))), : length of data.frame does not match number of Lines elementsDeserving
O
1

From SpatialPointsDataFrame to SpatialPolygonsDataFrame

library(sp)
library(raster)

### Example data: creating a SpatialPointsDataFrame object
x = c(1,2,5,4,3)
y = c(3,2,3,6,6)
df_points <- as.data.frame(cbind(x,y))
S <- SpatialPoints(cbind(x,y))
# S <- SpatialPoints(list(x,y))
# S <- SpatialPoints(data.frame(x,y))
S
plot(S)
spdf <- SpatialPointsDataFrame(S, df_points)
spdf
plot(spdf)
# crs(spdf) <- ("+proj=utm +zone=23 +south +datum=WGS84 +units=m +no_defs") ### add a crs

### Convert the SpatialPointsDataFrame to SpatialPolygons
(Sr1 = Polygon(spdf[,1:2]))
(Srs1 = Polygons(list(Sr1), "s1"))
(SpP = SpatialPolygons(list(Srs1), 1:1, proj4string= crs("+proj=utm +zone=23 +south +datum=WGS84 +units=m +no_defs"))) 
plot(SpP, col = 3:3, pbg="white", add=T) 
SpP ### can not write as shapefile

### Convert the SpatialPolygons to SpatialPolygonsDataFrame
shape_pol <- SpatialPolygonsDataFrame(SpP, match.ID=F, data= data.frame(x=spdf[1:1,1], y=spdf[1:1,2]))
shape_pol ### can be write as shapefile
plot(shape_pol, col = 4, add=T)

### write shapefile
library(rgdal)
writeOGR(shape_pol, paste0(getwd(), "/Output_shapes"), "p_to_shape_pol", driver="ESRI Shapefile")
Ober answered 10/3, 2021 at 13:44 Comment(0)

© 2022 - 2024 — McMap. All rights reserved.