3

I am new in analyzing raster datasets in R. I am trying to get Malawi's forest cover data at a lower administration level ( traditional authority level). So I downloaded shapefile at Traditional Authority level, and also got the Hansen data from Global Forest Change website (raster data).

I can extract the raster value at Malawi's border, but I am not getting at the traditional authority (polygon) level. This is my code, and I tried to explain what I did.

##Malawi Traditional Authority level Forest Cover Annual##

username <- "Malawi"
getwd()
setwd("")


library(sp)
library(rgdal)
library(raster)
library(gdalUtils)
library(parallel)
library(rgdal)


World= readOGR("C:/Users/howlade2/Documents/Malawi Maps/Malawi Forest/gadm28.shp", "gadm28")

head(world@data)

Malawi = world[world$NAME_0 == "Malawi",]

head(Malawi@data)

plot(Malawi)
##Raster Data from Hansen##

lossyear_malawi <- raster("C:/Users/howlade2/Documents/Malawi Maps/Malawi Forest/gadm28.shp/Hansen_GFC2015_lossyear_10S_030E.tif")


plot(lossyear_malawi, add=FALSE)
plot(Malawi, border = 'black', lwd = 4, add =TRUE)


##get values from raster##

time1 <- proc.time()
startTimeGB <-time1
clip1_gb <- crop(lossyear_malawi, extent (Malawi))
cat("crop1:","\n"); proc.time()-time1

time1 <- proc.time()
clip2_gb <- rasterize(Malawi, clip1_gb, mask =TRUE)
cat("crop2:","\n"); proc.time()-time1

time1 <- proc.time()
ext_gb <- getValues(clip2_gb)
cat("getValues:","\n"); proc.time()-time1

time1 <- proc.time()
tab_gb <- table(ext_gb)
cat("tabulate:","\n"); proc.time()-time1

time1 <- proc.time()
mat_malawi <- as.data.frame(tab_gb)
cat("as.data.frame(tab):","\n"); proc.time()-time1

timeGB <-proc.time() -startTimeGB
cat("Total time to extract raster cell values in Malawi:","\n"); proc.time()-time1
print(summary(timeGB))

length(ext_gb)

length(tab_gb)

GBlossyears <- mat_malawi$Freq[2:length(mat_malawi$Freq)]
GBlosscodes <- mat_malawi$Freq[2:length(mat_malawi$ext_gb)]

barplot(height =GBlossyears, names.arg = GBlosscodes, xlab ='Code for year of green loss', beside = TRUE )

##So this is giving me data for whole Malawi##

##But I need the data at the spatial polygon level (traditional authorities). And I couldn't manage what to do...##

rasterlist <- list("C:/Users/howlade2/Documents/Malawi Maps/Malawi Forest/gadm28.shp/Hansen_GFC2015_lossyear_10S_030E.tif") 
outlist <- list()
for (i in 1:length(rasterlist)) { #for each raster in rasterlist
  r <- raster(rasterlist[[i]]) #read element i of rasterlist into R
  val <- getValues(r) #get raster values
  m <- mean(val,na.rm=T) #remove NAs and compute mean
  outlist[[i]] <- c(rasterlist[[i]],m) #store raster path with mean
  return("complete")
}

This is the code:

# Query Malawi code from getData('ISO3') looking for the pattern "Malawi"
getData('ISO3')[grep(pattern = "Malawi", getData('ISO3')[,"NAME"]),]

malawiPolygon0 <- getData('GADM', country = 'MWI', level = 0)
malawiPolygon1 <- getData('GADM', country = 'MWI', level = 1)
malawiPolygon2 <- getData('GADM', country = 'MWI', level = 2)
malawiPolygon3 <- getData('GADM', country = 'MWI', level = 3)






# Download set of GFC tiles, there are two tiles. If you want to download both remove the subset option "[2]".
# (be careful, it's about 800 MB of data)
download_tiles(tiles, 
               output_folder = getwd(), 
               images = c("treecover2000", "loss", "gain", "lossyear", "datamask"), 
               data_year = 2015)

# Load rasters (only for tile 2)
treecover2000_malawi <- raster("Hansen_GFC2015_treecover2000_10S_030E.tif")
loss_malawi <- raster("Hansen_GFC2015_loss_10S_030E.tif")
gain_malawi <- raster("Hansen_GFC2015_gain_10S_030E.tif")
lossyear_malawi <- raster("Hansen_GFC2015_lossyear_10S_030E.tif")
datamask_malawi <- raster("Hansen_GFC2015_datamask_10S_030E.tif")

# Plot rasters Tiles for Malawi 
par(mfrow = c(3,2))

plot(treecover2000_malawi, main = "Global Forest Change - Tree cover", xlab = "Longitude", ylab = "Latitude")
plot(malawiPolygon2,  border = "#555753", col = NA, add = TRUE)

plot(loss_malawi, main = "Global Forest Change - Loss", xlab = "Longitude", ylab = "Latitude")
plot(malawiPolygon2,  border = "#555753", col = NA, add = TRUE)

plot(gain_malawi, main = "Global Forest Change - Gain", xlab = "Longitude", ylab = "Latitude")
plot(malawiPolygon2,  border = "#555753", col = NA, add = TRUE)

plot(lossyear_malawi, main = "Global Forest Change - Lossyear", xlab = "Longitude", ylab = "Latitude")
plot(malawiPolygon2,  border = "#555753", col = NA, add = TRUE)

plot(datamask_malawi, main = "Global Forest Change - Datamask", xlab = "Longitude", ylab = "Latitude")
plot(malawiPolygon2,  border = "#555753", col = NA, add = TRUE)

par(mfrow = c(1,1)) # default

# Assuming that the resolution level is for example "Admin level 1"
# Rasterize every polygon from Malawi administrative levels 1 SpatialPolygonsDataFrame object (malawiPolygon1)

# Crop GFC raster tile with country plygon (admin level = 2)
clip1_gb <- crop(lossyear_malawi, extent(malawiPolygon2))

# Rasterize and apply mask 
clip2_gb <- rasterize(malawiPolygon2, clip1_gb, mask = TRUE) # rasterize

# Get raster data for every polygon from admin level = 2
numberOfPolygons <- length(malawiPolygon2@polygons)
rasterList <- list()

for (i in 1:numberOfPolygons) {

  cat(paste("Cropping and rasterizing Polygon ", i, " - ", malawiPolygon2[i,]$NAME_2, "\n", sep = ""))
  r <- crop(clip2_gb, extent(malawiPolygon2[i,]))
  r <- rasterize(malawiPolygon2[i,], r, mask = TRUE)
  rasterList[[i]] <- r
  gc() # garbage collector to release memory

}

# Apply mean function to every raster belonging to every admin 1 polygon:
lossyearGFCAdmin1LevelMean <- lapply(rasterList, FUN = function(x) mean(values(x), na.rm = TRUE))

# Final Data frame
df <- data.frame("name" = malawiPolygon2$NAME_1, "lossyear_mean" = unlist(lossyearGFCAdmin1LevelMean))

> df
  • Do you want to extract forest cover per authority level which I suppose is a polygon? Then you need zonal statistics. In R use the zonal () function from the raster package. – Jens Mar 26 '17 at 12:49
  • What mean Traditional Authority level? What gadm level is that? 0, 1, 2, 3? – Guzmán Mar 27 '17 at 16:16
  • Yes, it is 2. (Name_2 in the attribute table) – Aparna Howlader Mar 28 '17 at 7:11
3

You can use the administration level you want to extract the data and make your analysis. In this example, I used Admin Level 2. Try the commented and reproducible example below:

Load libraries, get and plot data:

## Malawi Traditional Authority level Forest Cover Annual ##
# Load libraries
library('rgdal')
library('raster')
library('mapview')
library('gfcanalysis') # Tools for Working with Hansen et al. Global Forest Change Dataset

# Get Malawi SpatialPolygonsDataFrame for different administrative levels

# Query Malawi code from getData('ISO3') looking for the pattern "Malawi"
getData('ISO3')[grep(pattern = "Malawi", getData('ISO3')[,"NAME"]),]

malawiPolygon0 <- getData('GADM', country = 'MWI', level = 0)
malawiPolygon1 <- getData('GADM', country = 'MWI', level = 1)
malawiPolygon2 <- getData('GADM', country = 'MWI', level = 2)

# Plot Malawi
plot(malawiPolygon0, main = "Plot Malawi", xlab = "Longitude", ylab = "Latitude", col = "#000000", lwd = 2)
plot(malawiPolygon1, col = "#2E3436", lwd = 1.5, add = TRUE)
plot(malawiPolygon2, col = "#555753", lwd = 1, add = TRUE)
box()

PlotMalawi

# Example using the great mapview package
mapview(malawiPolygon0, layer = "Admin. level 0", color = "#000000") + 
  mapview(malawiPolygon1, layer = "Admin. level 1", color = "#729FCF") + 
  mapview(malawiPolygon2, layer = "Admin. level 2", color = "#AD7FA8") 

MalawiMapview

Get Hansen et al., Global Forest Cover data for year 2015:

# Calculate the GFC product tiles needed for a given AOI (Area Of Interest)
tiles <- calc_gfc_tiles(malawiPolygon0)

# Plot 
plot(tiles, col = "#EF2929", lwd = 2, main = "Plot Malawi + Tiles", xlab = "Longitude", ylab = "Latitude")
plot(malawiPolygon0,  col = "#555753", add = TRUE)
box()

Tiles

mapview(malawiPolygon0, layer = "Admin. level 0", color = "#000000") + mapview(tiles, layer = "Global Forest Cover Tiles for Malawi", color = c("#EF2929", "#EF2929")) # using mapview    

# Download set of GFC tiles, there are two tiles. If you want to download both remove the subset option "[2]".
# (be careful, it's about 800 MB of data)
download_tiles(tiles[2], 
               output_folder = getwd(), 
               images = c("treecover2000", "loss", "gain", "lossyear", "datamask"), 
               data_year = 2015)

# Load rasters (only for tile 2)
treecover2000_malawi <- raster("Hansen_GFC2015_treecover2000_10S_030E.tif")
loss_malawi <- raster("Hansen_GFC2015_loss_10S_030E.tif")
gain_malawi <- raster("Hansen_GFC2015_gain_10S_030E.tif")
lossyear_malawi <- raster("Hansen_GFC2015_lossyear_10S_030E.tif")
datamask_malawi <- raster("Hansen_GFC2015_datamask_10S_030E.tif")

# Plot rasters Tiles for Malawi 
par(mfrow = c(3,2))

plot(treecover2000_malawi, main = "Global Forest Change - Tree cover", xlab = "Longitude", ylab = "Latitude")
plot(malawiPolygon0,  border = "#555753", col = NA, add = TRUE)

plot(loss_malawi, main = "Global Forest Change - Loss", xlab = "Longitude", ylab = "Latitude")
plot(malawiPolygon0,  border = "#555753", col = NA, add = TRUE)

plot(gain_malawi, main = "Global Forest Change - Gain", xlab = "Longitude", ylab = "Latitude")
plot(malawiPolygon0,  border = "#555753", col = NA, add = TRUE)

plot(lossyear_malawi, main = "Global Forest Change - Lossyear", xlab = "Longitude", ylab = "Latitude")
plot(malawiPolygon0,  border = "#555753", col = NA, add = TRUE)

plot(datamask_malawi, main = "Global Forest Change - Datamask", xlab = "Longitude", ylab = "Latitude")
plot(malawiPolygon0,  border = "#555753", col = NA, add = TRUE)

par(mfrow = c(1,1)) # default

rasters

# Crop GFC raster tile with country polygon (admin level = 0)
clip1_gb <- crop(lossyear_malawi, extent(malawiPolygon0))

# Rasterize and apply mask 
clip2_gb <- rasterize(malawiPolygon0, clip1_gb, mask = TRUE) # rasterize

Calculate Lossyear mean for every Admin. Level 2 polygon:

# The resolution level is for "Admin level 2"
# Rasterize every polygon from Malawi administrative levels 2 SpatialPolygonsDataFrame object (malawiPolygon2)


# Get raster data for every polygon from admin level = 2 
numberOfPolygons <- length(malawiPolygon2@polygons)
rasterList <- list()

for (i in 1:numberOfPolygons) {

if (!is.null(intersect(extent(malawiPolygon2[i,]), clip2_gb))) {

cat(paste("Cropping and rasterizing Polygon ", i, " - ", malawiPolygon2[i,]$NAME_2, "\n", sep = ""))
r <- crop(clip2_gb, extent(malawiPolygon2[i,]))
r <- rasterize(malawiPolygon2[i,], r, mask = TRUE)
rasterList[[i]] <- r
gc() # garbage collector to release memory

} else {

cat(paste(" Attention: Polygon ", i, " - ", malawiPolygon2[i,]$NAME_2, " don't intersect with raster! \n", sep = ""))
rasterList[[i]] <- NA

}
}

# Apply mean function to every raster belonging to every admin 2 polygon:
lossyearGFCAdmin2LevelMean <- lapply(rasterList, FUN = function(x) mean(values(x), na.rm = TRUE))

# Final Data frame
df <- data.frame("name" = malawiPolygon1$NAME_2, "lossyear_mean" = unlist(lossyearGFCAdmin2LevelMean))

> df     

     name lossyear_mean
1    Balaka Town   0.010672820
2     TA Kalembo   0.017090464
3     TA Nsamala   0.023450509
4  Blantyre City   0.232012753
5     TA Chigaru   0.007986067
6    Balaka Town   0.065805384
7     TA Kalembo   0.054884431
8     TA Nsamala   0.025754032
9  Blantyre City   0.024921069
10    TA Chigaru   0.029292582

Note: I used only one tile, but you need the second tile for the north zone. If you don't use both tiles, some polygons won't intersect in the north zone of Malawi with the lossyear raster. You will need to merge rasters: merge(raster1, raster2)

# Crop GFC raster tile with country polygon (admin level = 0)
clip1_gb_tile1 <- crop(lossyear_malawi1, extent(malawiPolygon0)) # tile1
clip1_gb_tile2 <- crop(lossyear_malawi2, extent(malawiPolygon0)) # tile2

# Rasterize and apply mask 
clip2_gb_tile1 <- rasterize(malawiPolygon0, clip1_gb_tile1, mask = TRUE) 
clip2_gb_tile2 <- rasterize(malawiPolygon0, clip1_gb_tile2, mask = TRUE) 

clip2_gb_tile12 <- merge(clip2_gb_tile1, clip2_gb_tile1)

After doing that, you will use clip2_gb_tile12 as clip2_gb used in this example!

  • Thanks a lot! This is super helpful, I would not understand this much details! I followed all the steps by changing the admin level to TWO (traditional authority- there are 256). But this is giving error "Error in .local(x, y, ...) : extents do not overlap". I searched a bit- looks like projection might be wrong. Would you please help me a bit more? How to correct it ? – Aparna Howlader Mar 28 '17 at 7:00
  • @AparnaHowlader please see the answer again, I have updated the answer to 1) manage null intersections between rasters and polygons, 2) use admin level 2 and 3) give you some advices. Please consider vote up my question if it was useful for you and check as answered if it answers your question! =) – Guzmán Mar 28 '17 at 13:10
  • Yeah! That works perfect! THANKS SO MUCH, this post helped me to learn bunches of basic things! I really appreciate- – Aparna Howlader Mar 31 '17 at 5:20
  • One more question that I wanted to ask you- So here forestlossyear raster is actually giving us a annual indicator of forest loss. So if there is a forest loss in 2014, it has a value "14". If I want to extract that value per year (so I want to know the number of pixel that has a forest loss by year), I do not need the mean- I need something like count the number of pixel that has a value "14". I was trying to modify this line: lossyearGFCAdmin2LevelMean <- lapply(rasterList, FUN = function(x) mean(values(x), na.rm = TRUE)) But it looks like FUN doesn't have a count thing? – Aparna Howlader Mar 31 '17 at 5:22
  • For one unit, I get the thing I want from this: ext_gb <- getValues(rasterList[[2]]) I want to do this for all 256 units, and it might be close- but something is wrong... ext_gb_all <- matrix(,length(rasterList),ncol=14) for (i in 1:256) {ext_gb_all <- getValues(rasterList[i, 1:256]) } – Aparna Howlader Mar 31 '17 at 8:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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