For the work, we have to load the following libraries:

packages <- c("terra", "tidyverse", "dplyr")

Method

The process requires collecting district wise data on percentage of household having drinking water facilities and sanitation. The data is available from NHM data for 2015-2016 and 2019-2020. The data is saved in a csv format and merged to the district shapefile. Since the number of districts in Assam has changed over the years, the data for 2015-2016 is available for 27 districts while the data for 2019-2020 is available for 33 districts. The percentages are subtracted by 100 to get the percentage of households not enjoying drinking water and sanitation facilities. The calculated percentage is rasterised. this code is run exclusively for drinking water facilities. The same logic can be applied for sanitation facilities.

deprv_nhm_2015 <- read.csv("state/assam_district_27_nhm_drnkg_sanit.csv")
deprv_nhm_2019 <- read.csv("state/assam_district_33_nhm_drnkng_sanit.csv")

districts_27 <- vect("assam_district_27_4326.shp")
districts_33 <- vect("assam_district_33_4326.shp")

distr_nhm_27 <- merge(districts_27, deprv_nhm_2015, by='OBJECTID')
distr_nhm_33 <- merge(districts_33, deprv_nhm_2019, by='OBJECTID')

#reference image
var <- rast("state/pop_ind_ppp_2015_UNadj.tif")
plot(var)

#subtracting with 100 to get the deprived percentage
distr_nhm_27$drnkng_2015 <- (100-distr_nhm_27$drnkng_2015)/100
distr_nhm_33$drnkng_2019 <- (100-distr_nhm_33$drnkng_2019)/100

#rasterising the data
nhm_2015 <- rast(nrow = nrow(var), ncol = ncol(var))
ext(nhm_2015) <- ext(var)
res(nhm_2015) <- res(var)
crs(nhm_2015) <- crs(var)
nhm_2015 <- rasterize(distr_nhm_27, nhm_2015, 'drnkng_2015')
plot(nhm_2015)

nhm_2019 <- rast(nrow = nrow(var), ncol = ncol(var))
ext(nhm_2019) <- ext(var)
res(nhm_2019) <- res(var)
crs(nhm_2019) <- crs(var)
nhm_2019 <- rasterize(distr_nhm_33, nhm_2019, 'drnkng_2019')
plot(nhm_2019)

#user-defined functions to calculate maps of other years; assume the percentage to remain same for the gap years

norm_map_15 <- function(yr){
  var <- rast(sprintf("state/pop_ind_ppp_%s_UNadj.tif",yr))
  var <- var*nhm_2015
  return(var)
}

norm_map_19 <- function(yr){
  var <- rast(sprintf("state/pop_ind_ppp_%s_UNadj.tif",yr))
  var <- var*nhm_2019
  return(var)
}

#calling the user-defined functions
year <- 2015:2018
layermaps_15 <- sapply(year, norm_map_15)
year <- 2019:2022
layermaps_19 <- sapply(year, norm_map_19)

Finally, all the generated maps are saved.