[GIS] Need to extract rectangle coordinates from a shapefile

extractrshapefile

New R user. I'm trying to create an observation window for a point pattern in R using a shapefile imported from ArcGIS. I've imported the shapefile and it is recognized by R as a feature class type:polygon. I need the xmin, xmax, ymin, ymax of the extent and that is all in order to specify these in creating the owin object. I can't seem to extract these values from the polygon. I know this is very basic but I can't put 2 and 2 together here.

Best Answer

You can create a rectangle with the boundary box of your polygon. I use rgdal package to make a reproducible example, but you could use only raster package with shapefile() function instead of readOGR():

library(rgdal)
library(raster)

dsn <- system.file("vectors", package = "rgdal")[1]
ogrListLayers(dsn)
polygon <- readOGR(dsn=dsn, layer="scot_BNG")

Creates an extent object using bbox and convert it to polygon:

e <- extent(polygon@bbox)

e <- as(e,"SpatialPolygons")

See the results:

plot(polygon)
plot(e, add= T)

Related Solutions

[GIS] Reprojection of shapefile in R and final extraction of values

The above code looks more like you just clone shp into shputm and then assign the output of crs(shp) to shputm without performing an actual reprojection. Anyway, if you import both the shapefile and the NDVI raster, and then reproject shp using spTransform, subsequent data extraction should work out fine. Also, the output of extent(shp_utm) roughly agrees with the manually extracted shapefile extent you depicted above.

# Required packages
library(rgdal)
library(raster)

# Shapefile import
shp <- readOGR(dsn = "data", layer = "testflaechen_20131203")
# crs(shp)

# Raster import
rst <- raster("data/ndvi_LE71890272009095ASN00.tif")
# crs(rst)

# Shapefile reprojection
shp_utm <- spTransform(shp, crs(rst))

# compareCRS(rst, shp_utm)
# extent(shp_utm)

# Value extraction
extract(rst, shp_utm)

Process Vector to Raster Faster with R – Speed Optimization Techniques

I tried to "parallelize" the function rasterize using the R package parallel in this way:

split the SpatialPolygonsDataFrame object in n parts
rasterize every part separately
merge all the parts into one raster

In my computer, the parallelized rasterize function took 2.75 times less than the no-parallelized rasterize function.

Note: the code below download a polygon shapefile (~26.2 MB) from the web. You can use any SpatialPolygonDataFrame object. This is only an example.

Load libraries and example data:

# Load libraries
library('raster')
library('rgdal')

# Load a SpatialPolygonsDataFrame example
# Load Brazil administrative level 2 shapefile
BRA_adm2 <- raster::getData(country = "BRA", level = 2)

# Convert NAMES level 2 to factor 
BRA_adm2$NAME_2 <- as.factor(BRA_adm2$NAME_2)

# Plot BRA_adm2
plot(BRA_adm2)
box()

# Define RasterLayer object
r.raster <- raster()

# Define raster extent
extent(r.raster) <- extent(BRA_adm2)

# Define pixel size
res(r.raster) <- 0.1

Figure 1: Brazil SpatialPolygonsDataFrame plot

Simple thread example

# Simple thread -----------------------------------------------------------

# Rasterize
system.time(BRA_adm2.r <- rasterize(BRA_adm2, r.raster, 'NAME_2'))

Time in my laptop:

# Output:
# user  system elapsed 
# 23.883    0.010   23.891

Multithread thread example

# Multithread -------------------------------------------------------------

# Load 'parallel' package for support Parallel computation in R
library('parallel')

# Calculate the number of cores
no_cores <- detectCores() - 1

# Number of polygons features in SPDF
features <- 1:nrow(BRA_adm2[,])

# Split features in n parts
n <- 50
parts <- split(features, cut(features, n))

# Initiate cluster (after loading all the necessary object to R environment: BRA_adm2, parts, r.raster, n)
cl <- makeCluster(no_cores, type = "FORK")
print(cl)

# Parallelize rasterize function
system.time(rParts <- parLapply(cl = cl, X = 1:n, fun = function(x) rasterize(BRA_adm2[parts[[x]],], r.raster, 'NAME_2')))

# Finish
stopCluster(cl)

# Merge all raster parts
rMerge <- do.call(merge, rParts)

# Plot raster
plot(rMerge)

Figure 2: Brazil Raster plot