# How can I determine the area of each block of a SpatialGridDataFrame in R?

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
# ...
``````
• 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

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')
`````` 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) * res(r)
``````

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
``````

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