Unknown value while calculating land cover type buffer zones in ArcMap and r

78 views Asked by At

I am trying to calculate the percentages of land cover types within a 2.5km buffer zone through ArcMap and R. When I calculate the buffer zone TIFs in R using Landscapemetric, there is an unknown value that is calculated. The values are 10,30,40 and 50,60,80, which are all values present in the land cover data, but the last value is 128, which is not present. I believe it could be a background value, but I am not sure. If anyone knows how to remove this value so it does not show up on my calculations, that would be fantastic.

Land cover calculations buffer zone

I used landscapemetrics to calculate the land cover buffer zones which I thought would give the values shown on the buffer zones. The results gave me an extra unknown value.

2

There are 2 answers

0
flaviomcosta On

These values are appearing because the circle of the buffer will not be perfect due the square shape of the pixels. Rounded buffer boundaries tend to cut off edge pixels that turns it into NA's or other value.

Following an example of a cut raster by a buffer to show a missing boundaries pixels. enter image description here

To avoid these values you can calculate the percentages only for the valid classes that you are interested in your buffer.

Using only R you can extract the values of the raster layer with your buffer and count the number of pixels of each class in relation to total applying a simple rule of three to access the percentage of it. In this sense you calculate the percentage for all classes of which you are interest instead for the whole landscape (or buffer) thus you won't have a value that you are not interested.

I made an example with what you need to do in your dataset to achieve the goal of calculate landscape classes percentages within buffer cuts of a landscape raster image.

#packages
library(raster)
library(rgeos)

#an example of raster file with utm projection (to meters)
r = raster(nrow=10000, ncol=10000)
values(r) = sample(10, ncell(r), replace = TRUE)

#Important adjust the zone to your zone to project correctly
crs(r) = "+proj=utm +zone=20 +datum=WGS84 +units=m +no_defs"

#a point to create a buffer example 
pt <- SpatialPoints(data.frame(xx = 50, yy = 10))

#project the point to utm. Need to be in the same projection of the raster, utm to calculate the buffer in meters
proj4string(pt) = CRS("+proj=utm +zone=20 +datum=WGS84 +units=m +no_defs")

#making a buffer. In your case just read your buffer if already done (terra package is a good choice)
ptb  = gBuffer(pt, width = 25)

plot(r, legend = FALSE)
plot(ptb, add = TRUE)

enter image description here

#extract the raster values inside the buffer boundaries
buf_ex = unlist(extract(r, ptb))
#here we have only the valid classes values of the first buffer figure in this answer.

#now we construct a loop to calculate the percentages for each class (ps: maybe apply family offer a better option)
perc = data.frame(class = sort(unique(buf_ex)))
for(i in 1:nrow(perc)){
  #calculate the percentage of by the number of pixels
  perc[i, "percentage"] = length(buf_ex[buf_ex%in%i])*100/length(buf_ex)
}

The data.frame "perc" will show the percentages of each class within the buffer in the column "percentage".

If you are interested in the classes area or specific class area, simply multiply the pixel size by the number of pixels in the each class.

0
Robert Hijmans On

128 is probably the missing-value flag that is not identified as such. To exclude from it you could (A) register it, or (B) you could remove it altogether prior to your analysis.

library(terra)
#A
NAflag(x) <- 128
#B
y <- subst(x, 128, NA)

You could also simply remove the counts of 128 from your results.

You can extract and tabulate the values like this (with SpatVector pts)

b <- buffer(pts, 25000)
e <- extract(x, b, fun="table")

For example:

library(terra)
r <- rast(res=0.1)
set.seed(1)
values(r) <- sample(c(10,30,40,50,60,80), ncell(r), replace = TRUE)
pts <- vect(cbind(x=c(0, 20, 50), y=c(-20, 0, 10)), crs="+proj=longlat")

b <- buffer(pts, 25000)
e <- extract(r, b, "table")
e
#  ID 10 30 40 50 60 80
#1  1  2  1  4  2  5  2
#2  2  1  5  3  2  3  2
#3  3  3  3  2  5  0  3