Once the weights for the different classes of a layer are calculated using concepts borrowed from Bayesian statistics, we need to apply the weights to the original maps for further analysis. Her we show the codes for assigning the weights to indices - NDVI and NDBI, and land use map.
For the work, we have to load the following libraries:
packages <- c("terra","dplyr","tidyverse","caret")
The weights are evaluated for years from 2018 to 2022 when the raw data (Sentinel-2) is available is extract the indices and various classes of the land use map.
The steps for weighing the maps are as follows: (1) Read the table where the weights of each class of the variable are stored. (2) Read the map which is categorized into classes. For land use the original land use map is reclassified into 4 groups while for the others quantiles are used to classify the images into 5 groups. (3) Filter outs the weights of each variable from the table. (4) Reclassify the map in (2) by assigning the weights to the classes (5) Normalize the values in the map (6) For land use, the map in (5) is the final map which can be saved afer resampling and used for later analysis. (7) For the others, read the original map. (8) Take the product of map in (5) with map in (7). (9) Re-sample the result to the dimensions and extent of a selected map and save the results.
#RUN THE LOOP FOR ALL LAYERS
for (yr in 2018:2018){
#read the frequency ratio table
f <- read.table(sprintf('conditioningfactors/frequencyratio_Aug_%s',yr))
colnames(f) <- c('factor', 'Class', 'FR')
FR <- f
state_4326 <- vect("assam_State.shp")
FR$FR_norm <- FR$FR
#read the classified conditioning factors that varies across years
#LU - the original LU map with 11 classes is reclassified into 5 classes
lu <- rast(sprintf("conditioningfactors/classified_LU_Aug_%s.tif",yr))
plot(lu, main = "Landuse")
#NDBI AND NDVI - 2018 TO 2022
built <- rast(sprintf("conditioningfactors/indices_NDBI_Aug_%s.tif",yr))
plot(built, main = "Built density")
veg <- rast(sprintf("conditioningfactors/indices_NDVI_Aug_%s.tif",yr))
plot(veg, main = "Vegetation density")
#reclassify the classes with the FRs
FR_wt <- FR %>% filter(factor == "LU")
lu_FR <- classify(lu, cbind(id=FR_wt$Class, v=FR_wt$FR_norm))
lu_FR <- lu_FR/ max(lu_FR[], na.rm = T)
plot(lu_FR, main = "Land use - FR")
#reclassify the classes with the FRs
FR_wt <- FR %>% filter(factor == "NDBI")
built_FR <- classify(built, cbind(id=FR_wt$Class, v=FR_wt$FR_norm))
built_FR <- built_FR / max(built_FR[] , na.rm = T)
plot(built_FR, main = "Built density - FR")
#reclassify the classes with the FRs
FR_wt <- FR %>% filter(factor == "NDVI")
veg_FR <- classify(veg, cbind(id=FR_wt$Class, v=FR_wt$FR_norm))
veg_FR <- veg_FR / max(veg_FR[], na.rm = T)
plot(veg_FR, main = "Vegetation density - FR")
#as it is a categorical variable the FR is taken as the final value
lu_rl <- lu_FR
built_rl <- rast(sprintf("conditioningfactors/built_Actual_values_%s.tif",yr))
veg_rl <- rast(sprintf("conditioningfactors/vegetation_Actual_values_%s.tif",yr))
#applying weights to the layers
lu_rl_wt <- lu_FR
plot(lu_rl_wt, main = "WEIGHTED LANDUSE")
built_rl_wt <- built_rl * built_FR
plot(built_rl_wt, main = "Weighted built density")
veg_rl_wt <- veg_rl * veg_FR
plot(veg_rl_wt, main = "Weigted vegetation density")
slope <- rast("conditioningfactors/slope.tif")
dummy <- rast(ext=ext(slope),crs=crs(slope, proj=T, describe=FALSE, parse=FALSE),
nrow=dim(slope)[1],ncol=dim(slope)[2])
lu_rl_wt <- resample(lu_rl_wt, dummy)
built_rl_wt <- resample(built_rl_wt, dummy)
veg_rl_wt <- resample(veg_rl_wt, dummy)
#save the resampled ones
}
The same logic is followed for the demographic layers of socio-economic
variable namely, population, population density, distribution of
elderly, young population, sex ratio, AAY population, deprived
population, households with no phone of land line, households with no
drinking water or sanitation (code :
floodHaz_Method2_Aug_factors_changeacrossyears_for_demo_layer.R) and for
static layers of various flood conditioning factors (code :
frequencyratio_distancefromroad.R).