3

(I asked this question on stackoverflow, but since there were no solutions, and it has a GIS component, I thought I'd ask here).

I have an excel file with minimum and maximum projected x and y coordinates. I used the function rect to make rectangles

The code looks something like this, (the plot boundaries are different in my real script, but you get it hopefully):

op<- par(bg= "white")
plot(c(0,1000), c(0,1000), type = "n", xlab="", ylab= "", main= "poly test")
xmin=0;ymin=0;xmax=1000;ymax=1000
poly_test<- rect(xleft= xmin, ybottom=ymin, xright= xmax, ytop= ymax, lwd=3)

So far so good. But now I want to export it to QGIS or some other GIS application, or at the least have R recognize that these are spatial objects. When I try this,

poly_test2 <- SpatialPolygons(poly_test)

it gives an error. First that error was is.list(Srl) is not true. Then I turned poly-test2 into a list, and now I get "" Cannot get a slot('area") from an object of type "null"

What should I do? I've also tried setting the xmin, ymin etc as coordinates first, but it doesn't work either.

I also want to label my polygons based on an id listed in the excel file, how to do that in R? And will that label show up after exporting, in QGIS?

EDIT2: There was a projection in the projection file,.. I forgot to set something in QGIS , I think.. Anyway, I adapted the script by including a "for loop", since I have a spreadsheet with a bunch of minimum and maximun x and y coordinates. The for loop looks something like this:

for(i in 1: length(Excel_File$column))
{

e<- extent(Excel_File$XMIN[i],Excel_File$XMAX[i], Excel_File$YMIN[i],Excel_File$YMAX[i])
}

And then the rest of the script you suggested. When I apply plot(e), I see only one polygon. How do I get to see a plot of all of the polygons I created? Is my for loop not good? Do I need to crete e first with either e<-list(), or e<- data.frame etc?

  • I edit my answer to include your edit – aldo_tapia Sep 13 '17 at 18:34
3

First, rect() function is only graphical. It doesn't save any object.

Second, you need to create a SpatialPolygonsDataFrame and you'll be able to export polygons to QGIS and label it with Excel ID's. With the same data from your example:

# Load package
library(raster)
e <- extent(xmin,xmax,ymin, ymax) # extent object (a rectangle)

poly <- as(e,"SpatialPolygons")

Use proj4string() for CRS definition:

proj4string(poly) <- '+proj=utm +zone=21 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs '

poly <- SpatialPolygonsDataFrame(Sr = poly,data = data.frame(ID=1)) # data could be Excel column

Use rgdal package to export results:

library(rgdal)
# export it as shapefile
writeOGR(obj = poly,dsn = '/your/path', layer = 'layer_name', driver='ESRI Shapefile')

Process inside a for() loop

Your code:

for(i in 1: length(Excel_File$column))
{

e<- extent(Excel_File$XMIN[i],Excel_File$XMAX[i], Excel_File$YMIN[i],Excel_File$YMAX[i])
}

I suggest to use the following code:

library(raster)
library(rgdal)

proj <- '+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs ' # your CRS here

set.seed(111)

# Dummy file to make an reproducible example
Excel_File <- data.frame(column=1:4,
                         XMIN=sample(1:10,4),
                         XMAX=sample(1:10,4),
                         YMIN=sample(1:10,4),
                         YMAX=sample(1:10,4))

polys <- list() # create a list to store polygons

for(i in 1:length(Excel_File$column)){
  e <- extent(Excel_File$XMIN[i],Excel_File$XMAX[i], Excel_File$YMIN[i],Excel_File$YMAX[i])
  poly <- as(e,"SpatialPolygons")
  proj4string(poly) <-proj
  poly <- SpatialPolygonsDataFrame(Sr = poly,data = data.frame(ID=Excel_File$column[i]))
  polys[[i]] <- poly
}

# all polygons to an unique file

final <- do.call(bind,polys)

plot(final)

enter image description here

writeOGR(obj = final,dsn = '/your/path', layer = 'layer_name', driver='ESRI Shapefile')

# one file for each polygon

for(i in seq_along(polys)){ # you can include this in for loop above
  writeOGR(obj = polys[[i]],dsn = '/your/path', layer = paste0('layer_',i), driver='ESRI Shapefile')
}
  • I've tried to use this script... It does create the polygons etc.. However, when I add the layer to QGIS or ArcGIS, it doesn't recognize that I've already added a CRS system... I'm trying to use a projected coordinate system (the same projected coordinate system my coordinates are in), and I've read somewhere that this causes problems? – Minka Sep 11 '17 at 13:40
  • @Minka the output of writeOGR() are 4 files: name.shp, name.shx, name.dbf and name.prj. Check the last one with a text editor, it contains the CRS information in this format. If there is no .prj file, QGIS or ArcGIS will no recognise file's projection. – aldo_tapia Sep 11 '17 at 13:54
0

Your for loop will overwrite variable e; first make an empty vector or data.frame before your loop

e <- data.frame()

then, add a row to e each time the for loop iterates:

e <- rbind(e, extent(Excel_File$XMIN[i],Excel_File$XMAX[i], Excel_File$YMIN[i],Excel_File$YMAX[i])

then you make your rect()

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.