Basically, I computed a global distribution probability model in the form of ASCII, say:
gdpm
. gdpm
's values are all between 0 and 1.
Then I imported a local map from shape file:
shape <- file.choose()
map <- readOGR(shape, basename(file_path_sans_ext(shape)))
The next step, I rasterized gdpm
, and cropped using the local map:
ldpm <- mask(gdpm, map)
Then, I reclassified this continuous model into a discrete model (I divided the model into 6 levels):
recalc <- matrix(c(0, 0.05, 0, 0.05, 0.2, 1, 0.2, 0.4, 2, 0.4, 0.6, 3, 0.6, 0.8, 4, 0.8, 1, 5), ncol = 3, byrow = TRUE)
ldpmR <- reclassify(ldpm, recalc)
I've got a cropped and reclassified raster, now I need to summarize land cover, that is, to each level, I want to calculate its proportion of area in each region of the local map. (I don't know how to describe it in terminology). I found and followed an example (RobertH):
ext <- raster::extract(ldpmR, map)
tab <- sapply(ext, function(x) tabulate(x, 10))
tab <- tab / colSums(tab)
But I'm not sure if it works, since the output of tab
is confusing.
So how to compute land cover area correctly? How can I apply the correct method within each polygon?
My original data is too large, I can only provide an alternative raster (I think this example should apply a different reclassify matrix):
Example raster
Or you can generate a test raster (RobertH):
library(raster)
s <- stack(system.file("external/rlogo.grd", package="raster"))
writeRaster(s, file='testtif', format='GTiff', bylayer=T, overwrite=T)
f <- list.files(pattern="testtif_..tif")
I also have a question about plotting a raster:
r <- as(r, "SpatialPixelsDataFrame")
r <- as.data.frame(r)
colnames(r) <- c("value", "x", "y")
I do this conversion to make a raster plot-able with ggplot2, is there a more concise method?
The „area“ function of R's „raster package“ results in a vector of the cell sizes in the raster object (in km2), which differ from north to south. In order to calculate the entire area's size can be obtained by adding up all cell sizes or multiplying the median cell size by the number of cells (which I did here).
You can divide the area of the land use features by the area of the basin, then multiply by 100 to derive percent.
loki's answer is OK, but this can be done the raster way, which is safer. And it is important to consider whether the coordinates are angular (longitude/latitude) or planar (projected)
Example data
library(raster)
r <- raster(system.file("external/test.grd", package="raster"))
r <- r / 1000
recalc <- matrix(c(0, 0.05, 0, 0.05, 0.2, 1, 0.2, 0.4, 2, 0.4, 0.6, 3, 0.6, 0.8, 4, 0.8, 2, 5), ncol = 3, byrow = TRUE)
r2 <- reclassify(r, recalc)
Approach 1. Only for planar data
f <- freq(r2, useNA='no')
apc <- prod(res(r))
f <- cbind(f, area=f[,2] * apc)
f
# value count area
#[1,] 1 78 124800
#[2,] 2 1750 2800000
#[3,] 3 819 1310400
#[4,] 4 304 486400
#[5,] 5 152 243200
Approach 2. For angular data (but also works for planar data)
a <- area(r2)
z <- zonal(a, r2, 'sum')
z
# zone sum
#[1,] 1 124800
#[2,] 2 2800000
#[3,] 3 1310400
#[4,] 4 486400
#[5,] 5 243200
If you want to summarize by polygons, you can do something like this:
# example polygons
a <- rasterToPolygons(aggregate(r, 25))
Approach 1
# extract values (slow)
ext <- extract(r2, a)
# tabulate values for each polygon
tab <- sapply(ext, function(x) tabulate(x, 5))
# adjust for area (planar data only)
tab <- tab * prod(res(r2))
# check the results, by summing over the regions
rowSums(tab)
#[1] 124800 2800000 1310400 486400 243200
Approach 2
x <- rasterize(a, r2)
z <- crosstab(x, r2)
z <- cbind(z, area = z[,3] * prod(res(r2)))
Check results:
aggregate(z[, 'area', drop=F], z[,'Var2', drop=F], sum)
Var2 area
#1 1 124800
#2 2 2800000
#3 3 1310400
#4 4 486400
#5 5 243200
Note that if you are dealing with lon/lat data you cannot use prod(res(r)) to get the cell size. In that case you will need to use the area function and loop over classes, I think.
You also asked about plotting. There are many ways to plot a Raster* object. The basic ones are:
image(r2)
plot(r2)
spplot(r2)
library(rasterVis);
levelplot(r2)
More tricky approaches:
library(ggplot2) # using a rasterVis method
theme_set(theme_bw())
gplot(r2) + geom_tile(aes(fill = value)) +
facet_wrap(~ variable) +
scale_fill_gradient(low = 'white', high = 'blue') +
coord_equal()
library(leaflet)
leaflet() %>% addTiles() %>%
addRasterImage(r2, colors = "Spectral", opacity = 0.8)
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With