function to rasterize population on a grid

rasterize_FARM.R

# module load r-rgdal geos netcdf-fortran/4.5.2/gcc

# legge sezioni di censimento
library(sf)
st_read("/lustre/arpa/bonafeg/data/geo/Popolazione/FVGconSappada_SezCens.shp") -> ps

# legge popolazione
library(dplyr)
library(tidyr)
pp <- read.csv("/lustre/arpa/bonafeg/data/geo/Popolazione/Censimento2011/R06conSappada_indicatori_2011_sezioni.csv")%>%
 transmute(SEZ2011, population=P1)
left_join(ps,pp) %>% transmute(population=replace_na(population,0)) -> ps

# legge griglia
library(raster)
library("ncdf4")
raster("/lustre/arpa/bonafeg/scratch/bfm-source-apportionment/farm/base/FARM_conc_g1_20180101.nc")->gr
crs(gr) <- "+proj=utm +zone=33 +datum=WGS84 +units=km +no_defs "

# calcola frazioni di unita' censuaria su cella
library(exactextractr)
cc <- coverage_fraction(gr, ps, crop=F)

# rasterize
pop <- gr
values(pop) <- rep(0, ncell(gr))
pb <- txtProgressBar(min = 0, max = length(cc), style = 3)
for (i in 1:length(cc)) {
	pop <- pop + cc[[i]]*ps$population[i]/cellStats(cc[[i]],"sum")
	setTxtProgressBar(pb, i)
}
close(pb)
names(pop) <- "population"
crs(pop) <- "+proj=utm +zone=33 +datum=WGS84 +units=km +no_defs "

# plot
plot(pop)
cellStats(pop,"sum")

# salva
saveRDS(pop,"/lustre/arpa/bonafeg/data/geo/Popolazione/Pop_FARMFVG_2km.rds")