SpatialPolygons - Creating a set of polygons in R from coordinates

Question

I am trying to take create a set of polygons from vertex locations, saved in X,Y format.

Here is an example of my data - each row represents the vertices for one polygon. the polygons are squares

square <- rbind(c(255842.4, 4111578, 255862.4, 4111578, 255862.4, 4111558, 
                  255842.4, 4111558, 255842.4, 4111578, 255842.4, 4111578),
                c(257397.0, 4111309, 257417.0, 4111309, 257417.0, 4111289, 
                  257397.0, 4111289, 257397.0, 4111309, 257397.0, 4111309))

ID <- c("SJER1", "SJER2")'

I am using SpatialPolygons, thus my data need to be in a list. so i created a loop to attempt to get my data into a list format from a matrix.

I create a loop following code that I found in some other questions on this site. I broke out each step to try to understand why i am only getting one polygon as the output even through I have 2 sets of points.

for (i in 1:2) {  
  pts <- rbind(c(square[i,1], square[i,2]), c(square[i,3], square[i,4]), 
               c(square[i,5],square[i,6]), c(square[i,7],square[i,8]), 
               c(square[i,9],square[i,10]))
  sp1 <- list(Polygon(pts))
  sp2 <- list(Polygons(sp1,i))
  sp = SpatialPolygons(sp2)  
}
plot(sp)

Can you please help me understand how I adjust the code to write out two polygons instead of just one? And also, how I assign the ID to each polygon given I am using a matrix (square) as my starting dataset and if I assign a character id, it transforms all of my data into a character.

My ultimate goal is two polygons in the SpatialPolygons object, the first with the ID SJER1 and the second with the ID SJER2 stored in the SpatialPolygons object.

Then I will write this out to a shapefile.

Conceptualism

Answer 1

There's some information at ?'SpatialPolygons-class', but you more-or-less want to do the following:

polys <- SpatialPolygons(list(
  Polygons(list(Polygon(matrix(square[1, ], ncol=2, byrow=TRUE))), ID[1]),
  Polygons(list(Polygon(matrix(square[2, ], ncol=2, byrow=TRUE))), ID[2])
))

plot(polys)

The basic gist is that you need to create Polygon objects (e.g., from 2-column matrices with x coordinates in the first column and y coordinates in the second). These are combined in lists to create Polygons objects (each of which should have a unique ID). These Polygons objects are combined in a list to create a SpatialPolygons object. You can add a CRS if you like, with the proj4string argument to SpatialPolygons (see ?SpatialPolygons).

To write it out to an ESRI Shapefile, you'll need to convert it to a SpatialPolygonsDataFrame object by combining the polys object we created and some data. We'll just add the IDs as data for lack of anything more interesting.

polys.df <- SpatialPolygonsDataFrame(polys, data.frame(id=ID, row.names=ID))

and then write it out...

writeOGR(polys.df, '.', 'fancysquares', 'ESRI Shapefile')

The second argument ('.') says to write it out to the current working directory.

---

EDIT

To quickly create a SpatialPolygonsDataFrame when you have many rows describing polygons, you could use the following:

# Example data
square <- t(replicate(50, {
  o <- runif(2)
  c(o, o + c(0, 0.1), o + 0.1, o + c(0.1, 0), o)
}))
ID <- paste0('sq', seq_len(nrow(square)))

# Create SP
polys <- SpatialPolygons(mapply(function(poly, id) {
  xy <- matrix(poly, ncol=2, byrow=TRUE)
  Polygons(list(Polygon(xy)), ID=id)
}, split(square, row(square)), ID))

# Create SPDF
polys.df <- SpatialPolygonsDataFrame(polys, data.frame(id=ID, row.names=ID))

plot(polys.df, col=rainbow(50, alpha=0.5))

Trapezoid

thank you very much. i think i understand this better however, can you please tell me why either this: for (i in 1:2) { polys <- SpatialPolygons(list( Polygons(list(Polygon(matrix(square[i, ], ncol=2, byrow=TRUE))), ID[i]) )) } OR a <- vector('list', length(2)) for (i in 1:2) { a[[i]]<-list( Polygons(list(Polygon(matrix(square[i, ], ncol=2, byrow=TRUE))), ID[i]) ) } SpatialPolygons(a) don't work? i need to run this on many points so i will need some sort of loop. your time and explanation are much appreciated. i did review the documentation as well. – Conceptualism

@LWasser The first overwrites polys each time through the loop, so you end up with a SpatialPolygons object that only contains the last poly. The second will work if you remove the outermost list(), i.e. just do a[[i]] <- Polygons(list(Polygon(matrix(square[i, ], ncol=2, byrow=TRUE))), ID[i]). This is because a is already a list. – Trapezoid

@LWasser, see my edit for an approach that can be used to automate the process if you have many rows. (Also, apologies, I somehow forgot to include the code that creates the SPDF.) – Trapezoid

@jmaums thank you SO MUCH. I understand how this works now. The code works perfectly. I have a github page that i've been using to document these things as i figure them out / get help figuring them out or as i get code from others. i've thanked you on it as i sincerely appreciate your time in explaining / helping me -- Git hub URL to code – Conceptualism

@LWasser - you're welcome. Good luck with your research! – Trapezoid

it's really a shame sp doesn't have a function built in to handle this more succinctly – Austria