2

I'd like to be able to determine the area of each block inside a grid of (long, lat) coordinates. How can I do this? Alternatively, how can I determine the area of any Spatial Polygon with a long+lat projection? (I can probably extrapolate to my grid scenario if someone can show me how to do this.)

Setup

library(rgdal)
library(sp)

# US bounding box
usbb <- c(left=-125, bottom=20, right=-60, top=50)

# Generate the grid vertices (Each vertex is 10 degrees away from the one above it and the one to the left of it.)
usgridDF <- expand.grid(Lon=seq(usbb[["left"]], usbb[["right"]], by=10), Lat=seq(usbb[["bottom"]], usbb[["top"]], by=10))
usgridDF[, "GridId"] <- seq_len(nrow(usgridDF))

# Convert from data.frame to SpatialPointsDataFrame
usgrid <- SpatialPointsDataFrame(coords=usgridDF[, c("Lon", "Lat")], data=usgridDF[, "GridId", drop=F])

# Convert from SpatialPointsDataFrame to SpatialPixelsDataFrame
gridded(usgrid) <- TRUE

# Convert from SpatialPixelsDataFrame to SpatialGridDataFrame
usgrid <- as(usgrid, "SpatialGridDataFrame")
proj4string(usgrid) <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")

To Do

# Determine the area in each block of the grid, in sq km
# ...
  • 1
    areas computed directly from long-lat data can be obtained by geosphere::areaPolygon. All other methods need projected coordinates. – Edzer Pebesma Mar 1 '17 at 10:31
2

Quick solution is to convert to raster and use area() although I'm not 100% sure how accurate it'll be over big distances. Anyway:

library(raster)
usgrid_r <- raster(usgrid)
usgrid_r <- area(usgrid_r) # raster cell values are now area in sq km

# append as attribute to input grid
usgrid$area <- getValues(usgrid_r)

# this is just for plotting:
ugpolys <- as(usgrid, 'SpatialPolygonsDataFrame')
ugpolys$area_int <- as.integer(ugpolys$area)

library(tmap)
qtm(ugpolys, text = 'area_int')

plot of areas

2

Well, given your example you could coerce to a raster class object.

Add libraries and create data. Please note that I coerce the data to a SpatialPixelsDataFrame to skip some unnecessary code.

library(sp)
library(raster)

usbb <- c(left=-125, bottom=20, right=-60, top=50)
  usgridDF <- expand.grid(Lon=seq(usbb[["left"]], usbb[["right"]], by=10), Lat=seq(usbb[["bottom"]], usbb[["top"]], by=10))    
  usgridDF[, "GridId"] <- seq_len(nrow(usgridDF))

usgrid <- SpatialPixelsDataFrame(usgridDF[,c("Lon", "Lat")], data=usgridDF[, "GridId", drop=F])    
  proj4string(usgrid) <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")

We can then coerce the usgrid SpatialPixelsDataFrame to a raster class object. Since it is a grid (raster) than the area is uniform for each pixel. You can use the raster::res function to retrieve the x,y dimension of the pixels and get area. Just keep in mind that, since the data is not projected, the angular units are degrees. As such, you will need to convert degrees to an actual area. It is good practice to use projections that have associated planar coordinates (eg., meters, feet).

( r <- raster(usgrid, layer=1) )
res(r)[1] * res(r)[2]

You can also use the raster::area function to assign pixel areas to each pixel, as the raster value, that is computed using the latitudial and longitudinal span to approximate cell areas in square kilometers.

a <- area(r)
getValues(a)

You can, coerce the raster to a sp SpatialPolygonsDataFrame class object and pull the pixel areas from each polygon (just a good thing to know but completely unnecessary). This can also be done using rgeos::gArea, but not for non-projected geographic data.

r.poly <- rasterToPolygons(r, dissolve=TRUE )
sapply(slot(r.poly, "polygons"), function(x) sapply(slot(x, "Polygons"), slot, "area"))

r.poly@data$area <- getValues(a) #assign results from area function 
r.poly@data
0

Yet another alternative is to:

  1. convert to a SpatialPolygonDataFrame: as(., "SpatialPolygonsDataFrame")

  2. Project the data: spTransform(., CRS("+init=epsg:3410"))

  3. Use rgeos::gArea: gArea(., byid=TRUE)

Note: Contrary to the raster approach, you need to project the data, here I took a conical projection. Results differ slightly, but still have a correlation of 0.999

library(rgeos)

usgrid2 <- as(usgrid, "SpatialPolygonsDataFrame") 
usgrid3 <- spTransform(usgrid2, CRS("+init=epsg:3410"))

area <- gArea(usgrid3, byid=TRUE)

## compare (use data from previous post)
df <- data.frame(gArea=area/100000 , raster=getValues(usgrid_r))
cor(df)

Note: older versions of package sp do not allow for conversion , but there is the package Grid2Polygons that does it. Be careful that it does reorder the data according to id, so you will need to match it back.

library(Grid2Polygons)
usgrid2 <- Grid2Polygons(usgrid) 
area <- gArea(usgrid3, byid=TRUE)[usgrid$GridId]
  • good call on gArea, but check your version of sp, I had no trouble directly coercing SgridDF to SpolygonDF using as() with sp 1.2-4. – obrl_soil Mar 1 '17 at 11:10
  • oh, thanks for this pointing this out @obrl_soil ! Just edited accordingly. – Matifou Mar 1 '17 at 18:04

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.