Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Convert SpatialPointsDataFrame to SpatialLinesDataFrame in R

Tags:

r

spatial

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

like image 592
Richard Todd Avatar asked Jun 18 '14 11:06

Richard Todd


3 Answers

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
')
like image 57
Tomas Avatar answered Nov 17 '22 14:11

Tomas


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")
like image 34
Yuri Gelsleichter Avatar answered Nov 17 '22 14:11

Yuri Gelsleichter


## 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)
like image 4
mdsumner Avatar answered Nov 17 '22 15:11

mdsumner