Search code examples
rgisarcgisshapefile

Create rectangular polygon shapefiles with four corner coordinates

I have many observations and each of them has four corner coordinates to determine a rectangular area, such as x_min, x_max, y_min, and y_max. I want to create rectangular polygons using four corner coordinates and then export them as shapefiles for use in ArcGIS.

Here are some observations.

ID     xmax        xmin         ymax       ymin
1   -6987035.01 -8831307.63  6160489.27  4818866.55
2   -7457581.36 -7918649.51  5705536.08  5370130.4
3   -7527291.93 -7988360.08  5623289.83  5287884.15
4   -7632927.9  -7863461.98  5533399.89  5365697.05

Based on the coordinates, the left-lower point for the rectangle should be (xmin, ymin), the left-upper point for the rectangle should be (xmin, ymax), the right-lower point for the rectangle should be (xmax,ymin), and the right-upper point for rectangle should be (xmax, ymax).

How can I make a polygon file based on these four coordinates and export them as shapefiles which can be used in ArcGIS?

Solution

  • Try the following.

    I reformatted your data as a matrix of coordinates assuming that each row represents the bounds of a polygon.

    bounds <- matrix(
  c(
    -6987035.01, -8831307.63,  6160489.27,  4818866.55,
    -7457581.36, -7918649.51,  5705536.08,  5370130.4,
    -7527291.93, -7988360.08,  5623289.83,  5287884.15,
    -7632927.9,  -7863461.98,  5533399.89,  5365697.05
  ), nrow = 4, byrow = TRUE, dimnames = list(1:4, c('xmax','xmin','ymax','ymin'))
)

library(rgdal)
library(raster)

polygons_list = apply(bounds, 1, function(x){
  ## Loop throught rows to create extent objects
  out = extent(x[c(2, 1, 4, 3)]) ## x[c(2, 1, 4, 3)] to reoder the row to match the order required by raster::extent function. The order should xmin, xmax, ymin, ymax
  ## convert the extent objects to spatial polygons
  out = as(out, 'SpatialPolygons')
  ## convert the spatial polygons to spatial polygons dataframe
  out = as(out, 'SpatialPolygonsDataFrame')
  ## set CRS 
  crs(out) <- '+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs'
  out
})

## conbining the list of spatial polygons dataframes into one spatial polygons dataframe
polygons_list = do.call(rbind, polygons_list)
plot(polygons_list)

## export as shapefiles for future use in ArcGIS
writeOGR(polygons_list, getwd(), "polygons_list", driver="ESRI Shapefile", overwrite_layer=FALSE)