5

I have the following raster image that shows the different types of land covers on the island of Maui (there are 22 in total):

enter image description here

The island is divided into "mokus" or large districts, which can be seen in the following image:

enter image description here

How do I calculate how many pixels of each land cover type are within each moku from the land cover raster image above?

I've tried to do this by running the following code:

# Land cover raster image:
maui_lc <- raster::raster("hi_maui_2010_ccap_hr_land3.tif")

# Maui moku shapefile:
mokus <- rgdal::readOGR("Moku_Ridge_To_Reef_(DAR)")

maui_mokus <- crop(mokus, extent(-156.8, -155.9, 20.42, 21.116551948923718))

# Trying to extract raster pixel values within each moku:
test <- raster::extract(x = maui_lc, 
                        y = maui_mokus,
                        df = TRUE,
                        small = T,
                        method = "bilinear")

But my code just keeps running, so I don't know if this code actually produces the data I need. I've also gotten an error that says:

Error: cannot allocate vector of size 2.5 Gb

to which I've tried to solve by increasing my memory limit, but it still didn't work.

Improve this question
4
  • 2
    Forget using raster::extract function, it is notoriously slow and the package is slowly being replaced with terra. Instead, try the terra::extract function or even better exactextractr::exact_extract from the exactextractr package which can manage memory in a way that can process large problems. Just note that your polygon vector data must be in an sf class.
    – Jeffrey Evans
    Sep 24, 2021 at 4:20
  • What's the size of your raster in rows and columns? Have you tested your code on a small raster or with fewer districts? That's the best way of telling if your code is actually right to start with.
    – Spacedman
    Sep 24, 2021 at 6:22
  • What have you got method="bilinear" there? land cover classes shouldn't be interpolated (averaged) or you end up with classes like "3.25". Also, do you need small=TRUE? You'll get approximately the right answer without it for a big increase in speed.
    – Spacedman
    Sep 24, 2021 at 6:24
  • 1
    Same problem, different island: cran.r-project.org/web/packages/exactextractr/vignettes/…
    – dbaston
    Sep 26, 2021 at 19:31

3 Answers 3

Reset to default
1

I know this has been up for awhile and there area few comments but I don't see anything marked as resolved. Did you get something to work? If not, here is a way I came up with to get this done. It does use a for loop which is not in R best practices, but because Maui is relatively small and the LCLU raster I acquired was at 30-m spatial resolution, the code runs quickly and correctly.

rm(list=ls(all=T)) #clear memory
closeAllConnections() #close connections

library(rgdal)
library(raster)

#set working directory
setwd("your/working/directory")

#read in Moku shp file for Hawaii
#https://planning.hawaii.gov/gis/download-gis-data-expanded/
moku = readOGR(Sys.glob("moku.shp\\*.shp"))

#subset to get only the Moku in Maui
maui <- moku[moku$MOKUPUNI == "Maui",]

#read in the LCLU raster
#https://developers.google.com/earth-engine/datasets/catalog/USGS_GAP_HI_2001
#30-m raster clipped to Maui during gee export
lclu = raster("maui_lclu.tif")

#reproject shp file to match the raster 
maui = spTransform(maui,crs(lclu))

#create a blank dataframe to write into
jj = data.frame() 

for(i in 1:nrow(maui)){ #initate a loop for each moku in Maui
  tmp = maui[i,] #subset to get one moku
  crp = crop(lclu,tmp) #crop the lclu raster by moku
  msk = mask(crp,tmp) #mask out values outside of moku bounds
  tbl = data.frame(freq(msk)) #get the count of each lclu value
  tbl = tbl[!is.na(tbl$value), ] #remove NA values
  tbl$sqmeters = tbl$count*900 #add in a sq. meters column for extra info
  tbl$moku = tmp$MOKU #add in moku information as a column 
  jj = rbind(jj,tbl) #add rows to dataframe jj
  
  print(paste0("done with ",i," of ",nrow(maui)))} #give progress report for loop

jj #look at the dataframe

#export
write.csv(jj,file="lcluMauiMokuArea.csv")
Improve this answer
6
  • You can get the same thing using: 1) define classes classes = 1:10; 2) extract raster values and use lapply to get the counts v <- lapply(terra::extract(r, vect(p)), function(x) { table(factor(x, levels=classes)) } using the factor function allows us to set the levels which results in a consistent table when classes are absent in a polygon; 3) you can then use as.data.frame(do.call(rbind, v)) to get a data.frame with class pixel counts. You could actually do this all in one line.
    – Jeffrey Evans
    Feb 15, 2022 at 18:56
  • Great idea. Thank you for providing this solution. Seems you would need to know how many classes to declare...? For sure an improved solution to mine!
    – Zackary Werner
    Feb 15, 2022 at 19:21
  • You can get the unique values in your raster using terra::unique(x) which would be the "class" object in my above example. If you wrap prop.table around table it will return the proportions of each class.
    – Jeffrey Evans
    Feb 15, 2022 at 20:55
  • Okay interesting. I ran your suggestion and got some weird values considering there are 11 polygons in my shp file so 11 unique ids and land cover values range from 0 - 37 with only 25 of those classes present in the domain. My code creates four columns for the data frame, a value, count, area and polygon id. Your code supplies two rows, ID and landcover, for 25 columns which is the number of unique landcover types which makes sense but the values range from 0 - 250,000. Any idea what is causing this? I think maybe the counts aren't getting put into their own column...?
    – Zackary Werner
    Feb 15, 2022 at 21:49
  • 1
    I went ahead a posted an answer. The problem was terra now returning a data.frame and not a list so, tapply rather than lapply.
    – Jeffrey Evans
    Feb 16, 2022 at 0:07
1

Here is a solution using terra and sf (start migrating now, rgdal is gone starting Jan 01 2023 so, sp and raster days are numbered).

Add packages and create example data

library(terra)
library(sf)

lulc <- rast(xmin=571823.6, xmax=616763.6, ymin=4423540, 
          ymax=4453690, resolution=100, crs = "EPSG:26912")
  lulc[] <- sample(c(0:25), ncell(lulc), replace=TRUE)

s <- st_as_sf(spatSample(lulc, size=11, method="random", 
              as.points=TRUE))
  s <- st_buffer(s, 2500)             

plot(lulc)
  plot(st_geometry(s), add=TRUE)

We use terra::unique to pull the class values from our raster. Then, terra::extract to get the raster values associated with our polygons. Note that the object structure is a bit different that what raster returned and now we get a data.frame with an ID column, indicating the order of the polygon input and replicated to the number of pixels in the polygon (indexed using v$ID or v[,1]). The second (or more if multiband) column contains our raster values (indexed using v[,2]). So, tapply is now the suitable function to aggregate our results. The use of factor(x, levels=classes)) allows us to return a table of proportions representing all of the classes and not just the ones associated with a given observation. This syntax adds empty levels to the resulting factor vector so, all classes are returned in the results.

( classes <- sort(terra::unique(lulc)[,1]) )
  v <- terra::extract(lulc, vect(s))
    d <- as.data.frame(do.call(rbind,tapply(v[,2], v$ID, function(x) { 
           prop.table(table(factor(x, levels=classes)))})))
  names(d) <- paste0("class", names(d))

We can now pull some information from the class proportions like, which class is the majority and what is its proportion and assign them to our polygons.

s$class.maj <- apply(d, 1, function(x) names(d)[which.max(x)])
s$class.prop <- apply(d, 1, function(x) x[which.max(x)] )

plot(s["class.maj"])

However, my go to for polygon/raster extraction is still exactextracr::exact_extract because it returns fractional intersection of each pixel (allowing for weighted sum and mean). This function still returns a list object so, one would use lapply or for as your iteration function.

Improve this answer
0

This has been simplified in terra version 1.7.21. That is currently the development version. You should be able to install that version with

install.packages('terra', repos='https://rspatial.r-universe.dev')

Example data

library(terra)
#terra 1.7.21
v <- vect(system.file("ex/lux.shp", package="terra"))
r <- rast(system.file("ex/elev.tif", package="terra"))
r <- round((r-50)/100)
levels(r) <- data.frame(id=1:5, name=c("forest", "water", "urban", "crops", "grass"))

You can use extract to get cell counts (use exact=TRUE to consider cell fractions).

e <- extract(r, v, fun="table", na.rm=TRUE, exact=FALSE)

# add the regions label
data.frame(NAME_2=v$NAME_2[e[,1]], e[,-1])
#             NAME_2 forest water urban crops grass
#1          Clervaux      0     0    28   459    74
#2          Diekirch      2   123   200    66     3
#3           Redange      0   109   170   170    17
#4           Vianden      0    31    39    57     3
#5             Wiltz      0     1   161   303     8
#6        Echternach     14    76   233     1     0
#7            Remich     50   147    24     0     0
#8      Grevenmacher     19   221   137     2     0
#9          Capellen      0    25   305     0     0
#10 Esch-sur-Alzette      0   190   229    15     0
#11       Luxembourg      0   184   223    16     0
#12           Mersch      0   167   239    14     0

To compute areas, you can use expanse or zonal. First with expanse:

# rasterize the zones
zone <- rasterize(v, r, "NAME_2")
expanse(r, unit="km", byValue=TRUE, zones=zone, wide=TRUE)[,-1]

#               zone    forest       water     urban      crops     grass
#1          Diekirch  1.110530  68.3111509 111.08929  36.621714  1.662192
#2        Echternach  7.784033  42.2720014 129.61202   0.555644  0.000000
#3      Grevenmacher 10.591752 123.1872271  76.33934   1.114409  0.000000
#4            Remich 27.951955  82.1601534  13.41321   0.000000  0.000000
#5          Capellen  0.000000  13.9496106 170.14699   0.000000  0.000000
#8  Esch-sur-Alzette  0.000000 106.2479147 128.04839   8.392079  0.000000
#10       Luxembourg  0.000000 102.7248564 124.44157   8.921887  0.000000
#11           Mersch  0.000000  92.9318291 133.02760   7.792452  0.000000
#12          Redange  0.000000  60.6687117  94.55675  94.459184  9.444716
#14          Vianden  0.000000  17.1997610  21.63354  31.602716  1.662192
#15            Wiltz  0.000000   0.5546965  89.26629 168.018372  4.437191
#17         Clervaux  0.000000   0.0000000  15.49914 253.829542 40.900555

Now with zonal

x <- cellSize(r, unit="km")
zonal(x, c(r, zone), fun="sum", wide=TRUE)

#              group       water     urban      crops     grass    forest
#1          Capellen  13.9496106 170.14699   0.000000  0.000000  0.000000
#3          Clervaux   0.0000000  15.49914 253.829542 40.900555  0.000000
#6          Diekirch  68.3111509 111.08929  36.621714  1.662192  1.110530
#11       Echternach  42.2720014 129.61202   0.555644  0.000000  7.784033
#15 Esch-sur-Alzette 106.2479147 128.04839   8.392079  0.000000  0.000000
#18     Grevenmacher 123.1872271  76.33934   1.114409  0.000000 10.591752
#22       Luxembourg 102.7248564 124.44157   8.921887  0.000000  0.000000
#25           Mersch  92.9318291 133.02760   7.792452  0.000000  0.000000
#28          Redange  60.6687117  94.55675  94.459184  9.444716  0.000000
#32           Remich  82.1601534  13.41321   0.000000  0.000000 27.951955
#35          Vianden  17.1997610  21.63354  31.602716  1.662192  0.000000
#39            Wiltz   0.5546965  89.26629 168.018372  4.437191  0.000000

To get the ("exact") area with extract, you could do

e <- extract(c(x,r), v, exact=TRUE, na.rm=TRUE) |> na.omit()
e$area <- e$area * e$fraction
a <- aggregate(e[, "area", drop=FALSE], e[, c("ID", "name")], sum)
a <- data.frame(NAME_2 = v$NAME_2[a[,1]], a)[,-2]

reshape(a, idvar="NAME_2", timevar="name", direction="wide")
#             NAME_2 area.forest  area.water area.urban  area.crops area.grass
#1          Diekirch    1.110530  67.6791652  109.78461  37.3792518   1.835912
#2        Echternach    8.224001  42.6656508  130.98765   0.3979812         NA
#3            Remich   26.893107  81.6930623   13.35597          NA         NA
#4      Grevenmacher    9.685963 122.5968312   74.52150   1.1144093         NA
#6           Redange          NA  60.6745706   93.69050  93.9922323   9.444456
#7           Vianden          NA  17.0706616   22.18677  30.8643367   1.454355
#8             Wiltz          NA   0.7361318   89.07461 167.5479682   4.768619
#12         Capellen          NA  14.7992764  169.55486          NA         NA
#13 Esch-sur-Alzette          NA 105.7853834  127.80048   8.2993843         NA
#14       Luxembourg          NA 103.0913118  125.09963   8.9218868         NA
#15           Mersch          NA  92.6554656  132.88209   7.7924516         NA
#16         Clervaux          NA          NA   13.88577 251.4770800  40.600391
Improve this answer

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.