Summarizing values across polygons, for each layer in Raster Stack - How to merge values with shapefile?

Question

I have a raster stack with 26 layers that represent different weeks of the year. Within this area there are 854 irregularly shaped polygons that I am interested in summary values for (mean and sd), for each of the layers.

I believe that I have successfully extracted these values. However, when I try to merge them back with the polygon file, I am unsuccessful. I believe that the extract() function has stripped the values of their identifying criteria (Name_2). I tried converting the "Large list" output from the extract() function into a data frame, but that didn't work either.

Any suggestions?

Thank you.

Brazil0 <- getData('GADM', country='BRA', level=0)
Brazil1 <- getData('GADM', country='BRA', level=1)
Brazil2 <- getData('GADM', country='BRA', level=2)

MinasGerais<-subset(Brazil2, NAME_1 =='Minas Gerais')
e <- extent(MinasGerais) 
MG <- as(Brazil2[Brazil2$NAME_1 == 'Minas Gerais',], 'SpatialPolygons') 
row.names(MG) <- as.character(1:length(MG))

setwd("H:/Brazilian Arbovirus Project/ratsers_EpiWeek")
Precip_2017 <- stack("Precip_2017.tif")

Precip_crop <- crop(Precip_2017, e, snap="out")  
crop<- setValues(Precip_crop, NA) 
Precip_raster <- rasterize(MinasGerais, crop)   
Precip_MG<- mask(x=Precip_crop, mask=Precip_raster)

MG_Municip_mean_2 <- extract(Precip_MG, 
                           MG,
                           method='simple',
                           match.ID=FALSE, 
                           FUN=mean,
                           #sp=TRUE,
                           small=TRUE)  

output=data.frame(MG_Municip_mean_2)

MG_1 <- merge( MinasGerais, MG_Municip_mean_2, by='NAME_2', all=TRUE)
writeOGR(MG_1, getwd(), "MG_Municipalities", driver="ESRI Shapefile", 
check_exists=TRUE, 
overwrite_layer=TRUE)

When you say you "believe you have extracted these values", have you confirmed this by debugging the code, or are you just assuming it works? — FluffyKitten
can you add some example data using dput (of Precip_2017) ore create similar data? which packages are you using? — minem

aldo_tapia aldo_tapia · Accepted Answer · 2017-08-16T14:45:34

For this kind of calculation I prefer to rasterize polygons rather than use extract(), it's really easy to use this approach for multilayer rasters.

With a reproducible example:

library(sp)
library(raster)
library(dplyr)

Brazil0 <- getData('GADM', country='BRA', level=0)
Brazil1 <- getData('GADM', country='BRA', level=1)
Brazil2 <- getData('GADM', country='BRA', level=2)

set.seed(123)

r <- raster(ext = extent(Brazil0@bbox))

rList <- list()

for(i in 1:26){
  rList[[i]] <- setValues(r,values = rnorm(n = ncell(r)))
}

rs <- stack(rList)

The code above create a multilayer raster with random values, now it's time to rasterize and create a data.frame to made calculations:

B0r <- rasterize(Brazil0,r,'OBJECTID')
B1r <- rasterize(Brazil1,r,'OBJECTID')
B2r <- rasterize(Brazil2,r,'OBJECTID')

df <- as.data.frame(stack(B0r,B1r,B2r,rs))

names(df)[1:3] <- paste0('level',0:2)

names(df)[4:29] <- paste0('week',1:26)

I'll use Brazil2 polygons for this example, in this case there are polygons smallest than 1 pixel, so not all polygons will have information (maybe you're working with a better spatial resolution). Using dplyr package, you can compute easily any kind of formula for each polygon and merge back to the attributes table:

# example for level2

Brazil2@data <- df %>% select(-c(level0,level1)) %>% group_by(level2) %>%
  summarise_each(funs(mean(., na.rm = TRUE),sd(., na.rm = TRUE))) %>%
  left_join(x=Brazil2@data,by=c('OBJECTID'='level2'))

head(Brazil2)
##   OBJECTID ID_0 ISO NAME_0 ID_1 NAME_1 ID_2          NAME_2 HASC_2 CCN_2 CCA_2    TYPE_2    ENGTYPE_2
## 1        1   33 BRA Brazil    1   Acre    1      Acrelândia           NA       Município Municipality
## 2        2   33 BRA Brazil    1   Acre    2    Assis Brazil           NA       Município Municipality
## 3        3   33 BRA Brazil    1   Acre    3       Brasiléia           NA       Município Municipality
## 4        4   33 BRA Brazil    1   Acre    4          Bujari           NA       Município Municipality
## 5        5   33 BRA Brazil    1   Acre    5        Capixaba           NA       Município Municipality
## 6        6   33 BRA Brazil    1   Acre    6 Cruzeiro do Sul           NA       Município Municipality
##   NL_NAME_2 VARNAME_2  week1_mean  week2_mean week3_mean  week4_mean week5_mean  week6_mean  week7_mean
## 1                      0.02799743  0.65684118 0.18032148  0.15386044  0.5635764  0.13131793  0.66844905
## 2                      0.62831026 -0.46710306 0.56323394  0.16482915 -0.1190538 -0.42822608  0.11841901
## 3                     -0.04441109  0.30017178 0.24697659 -0.44659707  0.1560986  0.08040939  0.11568072
## 4                      0.04773195 -0.18690385 0.03840785 -0.06729135  0.3124202 -0.04866495 -0.38620625
## 5                     -0.18407894  0.21761910 0.30628065  0.11929847 -0.5596950 -0.45436160  0.34162667
## 6                     -0.31508001 -0.07326906 0.19979837 -0.07487962  0.2475976  0.13038549  0.08667631
##    week8_mean  week9_mean week10_mean week11_mean week12_mean week13_mean week14_mean week15_mean
## 1  0.29398763 -0.42470630 -0.39871869 -0.11133687 -0.11600877  0.04953248   0.2538490  -0.6387728
## 2  0.43123080 -0.19684622 -0.52333000  0.48386358  0.30651641  0.02512788  -0.0629341  -0.3729230
## 3 -0.36615371  0.26210474 -0.30212885 -0.15083871 -0.07095279 -0.11673077  -0.4448210   0.2400551
## 4 -0.12058451 -0.07914649  0.07276792  0.08345732 -0.17132937 -0.25823307  -0.1630691  -0.3555612
## 5 -1.28023743 -0.43201474 -0.92807227  0.32306609  0.76073801  0.50142559   0.7291509   0.6027042
## 6  0.08928839 -0.23224347  0.10913390  0.35135734  0.04670324  0.09490760   0.3855241   0.1820617
##   week16_mean week17_mean week18_mean  week19_mean week20_mean week21_mean week22_mean week23_mean
## 1  0.62806442 -0.23251001  0.41445121 -0.357856979  0.06260211  0.25732954 -0.08781137  0.33854357
## 2 -0.45495979 -0.05871667 -0.09732541 -0.131669502 -0.08885639 -0.48216217 -0.11940655 -0.52798951
## 3  0.31939261 -0.19266620  0.08671799  0.241000093  0.14864434  0.07044010  0.47460312  0.07706955
## 4  0.08420384 -0.35656266 -0.18245677  0.008064379  0.07319371  0.09271722 -0.13242697  0.29805436
## 5  0.62892059 -0.07006391  0.41960685  0.167582501  0.10296624  0.07797158  0.03335396  0.20283157
## 6  0.26078806  0.16028080  0.31722646 -0.166467152  0.25832317 -0.13020199 -0.05083594 -0.07325540
##   week24_mean week25_mean week26_mean  week1_sd  week2_sd  week3_sd  week4_sd  week5_sd  week6_sd
## 1   1.0375927  0.05829013 -0.08900512 0.5039079 0.7364619 0.6436874 1.1242115 0.8373088 1.2155468
## 2   0.5672509  0.44157869 -0.45395791 0.6429615 0.6969150 1.2544599 1.1263861 0.8122241 0.9647385
## 3  -0.3314678 -0.16537211  0.27106228 1.2177288 0.7593655 0.9915033 1.0651381 0.8324841 1.1640994
## 4   0.0992971  0.46187201  0.19528712 1.0670819 0.5727231 0.9797493 0.9037118 1.1892758 0.6578113
## 5  -0.5825252 -0.24512778  0.48153479 1.3876316 1.0093498 0.6523201 0.6672334 0.2975399 0.6660514
## 6  -0.1356949 -0.11380037 -0.07519057 0.8528147 1.0521164 0.6692769 0.9140993 1.1120830 0.9439429
##    week7_sd  week8_sd  week9_sd week10_sd week11_sd week12_sd week13_sd week14_sd week15_sd week16_sd
## 1 1.0242434 1.0283737 0.7440292 0.7977511 0.5330551 1.1060657 1.0586667 0.9557767 1.2021589 0.8015805
## 2 0.8650897 1.1318929 0.7925099 0.8416769 0.6735254 1.2130762 1.1684757 1.4122850 1.1372183 1.2514929
## 3 1.0337027 0.9709125 0.7129832 0.7574105 1.3073536 0.8376609 1.1161002 0.8083519 1.2163283 0.9789011
## 4 0.7550593 0.6019020 1.0640847 1.1187340 0.9331608 0.9246086 1.0412971 0.9577963 1.0302680 0.8863142
## 5 1.2286351 0.8620402 0.5704018 1.5840274 0.8929725 0.6766018 0.4504001 0.4747059 0.6959029 1.8061449
## 6 1.1981396 1.0422314 1.0996494 0.7749073 0.7755682 1.2341215 0.9699409 0.8534397 0.8832001 1.1267263
##   week17_sd week18_sd week19_sd week20_sd week21_sd week22_sd week23_sd week24_sd week25_sd week26_sd
## 1 1.7168281 1.3058866 0.7613507 1.4192652 1.0980767 0.5848721 1.1768701 1.2879231 1.1683099 0.8173507
## 2 0.8853379 1.2223339 1.0174595 0.8165643 0.7453179 1.4108344 1.4206444 1.0326178 0.5363038 1.1400113
## 3 0.7980645 1.2271442 1.1425325 0.9385259 1.1597657 1.0255039 1.1598827 0.7623564 0.7854232 0.8780667
## 4 0.9595817 1.1212473 0.8296976 1.1852795 1.2297023 1.0755448 1.0810457 0.8593361 1.0980879 0.7695382
## 5 1.3182950 0.4643851 0.8796118 1.0272719 0.6404252 0.4710161 0.5726681 0.6737512 1.3994816 1.0740231
## 6 1.1344667 0.9216528 1.3398224 0.7941615 1.2181886 1.1101733 1.0151286 1.3209190 1.0377582 1.0426156

Summarizing values across polygons, for each layer in Raster Stack - How to merge values with shapefile?

1 Answers