CODEWITHSUNDEEP
rextractrasterterra
Alex White
Alex White

Reputation: 31

Calculate polygon area for each land cover type (using spatraster object)

I'm trying to calculate the area (and ideally percentage) of a 6 bands within a spatraster (land.sub). I made this spatraster by cropping an original land cover raster to polygon data (which represent presence of seagrass). I'm trying to work out how much of the 6 land cover types are present on my spat raster (land.sub) and have tried using exact_extractr and played around with terra's extract function but can't seem to work out how to go about it. This is my code:

library(terra)
library(sf)
library(dplyr)
library(rgdal)
library(ggplot2)
library(exactextractr)
##READ IN RASTER###
rcp19<-rast("7landtypes/SSP1_RCP19/global_SSP1_RCP19_2050.tif")
crs(rcp19, proj = T)

###READ IN POLYGON DATA AND TRANSFORM TO SAME CRS AS RASTER(CEA) BY TYPING >CMU AND COPY PASTING CRS DETAILS IN OUTPUT#### 
cmu <- st_read("units-attributes-wgs84L2.gpkg") %>% 
  st_transform("+proj=cea +lat_ts=30 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs") %>% 
  
  ###FILTER FOR JUST AUSTRALIA FROM TERRITORY1 VARIABLE#####
filter(TERRITORY1 == 'Australia')

###FILTER TO JUST SEAGRASS###
sg<-cmu[, names(cmu)%in% c ("seagrass", "unit_ID")] %>% filter(seagrass==1)
plot(sg)

rcp19
######################EXACT EXTRACTR############################################

###CREAT OBJECT CALLED LAND.SUB WHICH COMBINES RCP19 ONLY WITH WHERE SEAGRASS PRESENT###
land.sub <- terra::crop(rcp19, vect(sg), mask = T)
plot(land.sub)

####RENAME NUMERIC VARIABLES TO LANDCOVER TYPES####
levels(land.sub)<-data.frame(id=1:6, cover = c("water", "forest", "grassland", "barren", "cropland", "urban"))
plot(land.sub)

I should also mention when attempting calculations with exact_extractr i kept getting error: columns undefined. Image of my map attached.

landsub raster with land cover raster and seagrass polygons cropped together

Upvotes: 1

Views: 1191

Answers (1)

Robert Hijmans
Robert Hijmans

Reputation: 47071

I cannot follow your example, but if you want to compute the area for different classes you can do that like I show below.

Example data

library(terra)
r <- rast(system.file("ex/elev.tif", package="terra"))
r <- round(r/100)
v <- vect(system.file("ex/lux.shp", package="terra"))

If you want to subset an area using polygons you can do something like

r <- crop(r, v[1:2, ], mask=TRUE)

You can count the number of cells in each class with freq

freq(r)
#  layer value count
#1     1     2    47
#2     1     3   219
#3     1     4   314
#4     1     5   475

To compute the area of each class:

a <- cellSize(r, unit="km")
zonal(a, r, fun=sum)
#  elevation     area
#1         2  26.1040
#2         3 121.6354
#3         4 174.0550
#4         5 262.7002

Upvotes: 3

Related Questions