data-raw/23-elabora_confini.R

In master-info/masteRomi: Dati, Poligoni e altro per analisi di Zone OMI

#############################################
# ELABORA CONFINI *ULTIMO PERIODO* ZONE OMI #
#############################################

masteRfun::load_pkgs('data.table', 'sf', 'rmapshaper')
fn <- max(list.files(file.path(bnd_path, 'OMI'), '20[0-9]'))

message('Leggo confini OMI ultimo periodo in WGS84')
yb <- qs::qread(file.path(bnd_path, 'OMI', fn))

message('Carico tabelle ancillari ')
yo <- fread('./data-raw/csv/zone.csv', select = c('catasto', 'termine', 'OMI'))
yo <- masteRgeo::comuni_storico[, .(CMN, catasto)][yo[is.na(termine)], on = 'catasto'][, termine := NULL]
yo <- masteRgeo::comuni[, .(CMN, RGNd)][yo, on = 'CMN']
setorderv(yo, c('RGNd', 'CMN'))

message('Applico "variazioni" a `zone` e `confini`')
# check `yo[!OMI %chin% yb$OMI]` vs `yb |> subset(!OMI %in% yo$OMI)`

message('Salvo versione semplificata nel pacchetto ')
k <- 0.2
ys <- unique(yo$RGNd) |> 
        sort() |> 
        lapply(\(x) yb |> subset(OMI %in% yo[RGNd == x, OMI]) |> ms_simplify(k, keep_shapes = TRUE, sys = TRUE)) |> 
        masteRgis::st_unisci_poli() |> 
        rbind(yb |> subset(!OMI %in% yo$OMI)) |>
        st_transform(4326)
if(nrow(ys) != nrow(yb)) stop('Alcune zone sono state perse nel processo di semplificazione. Aumentare il parametro `k`.')
st_write(ys, file.path(bnd_path, '20.gpkg'), 'OMI', append = FALSE, quiet = TRUE)
salva_dts_pkg(ys, 'OMI', as_db = FALSE, csv_in_pkg = FALSE)
if(file.info('./data/OMI.rda')$size > 1e5*1024)  stop('Il file `OMI.rda` supera 100MB. Diminuire il parametro `k`.')

message('Dissolvo i confini <OMI> per creare i confini Comunali <CMN>')
ysc <- ys |> merge(yo[, .(CMN, OMI)]) |> masteRgis::stm_dissolvi('CMN')
salva_dts_pkg(ysc, 'CMNo', as_db = FALSE, csv_in_pkg = FALSE)
st_write(ysc, file.path(bnd_path, '20.gpkg'), 'CMNo', append = FALSE, quiet = TRUE)

message('Inizio la procedura di associazione SZN>OMI.')

message(' - leggo confini sezioni con proiezione... ')
yz <- st_read(file.path(bnd_path, '0p.gpkg'), 'SZN', quiet = TRUE) |> subset(select = SZNid)
yb <- st_read(file.path(bnd_path, '0p.gpkg'), 'OMI', quiet = TRUE)

message(' - determino sezioni interamente coperte... ')
yc <- yz |> st_covered_by(yb)
for(x in which(sapply(yc, length) > 1)) yc[[x]] <- yc[[x]][1]
yc <- as.numeric(yc)
yc <- data.table( yz |> st_drop_geometry(),  yb[yc, 'OMI'] |> st_drop_geometry() )

message(' - determino sezioni intersecanti... ')
yc1 <- yc[!is.na(OMI)][, `:=`( area_int = 0, copertura = 100, N = 1 )]
yc2 <- yc[is.na(OMI)]
yc2 <- st_intersection(yz |> subset(SZNid %in% yc2$SZNid), yb)

message(' - determino superficie di intersezione... ')
yc2t <- data.table( 
            yc2 |> st_drop_geometry(), 
            area_int = as.numeric(st_area(yc2)) 
        )[, area_szn := sum(area_int), SZNid][, copertura := round(area_int / area_szn * 100)]
yc2t <- yc2t[order(SZNid, -copertura)][, N := 1:.N, SZNid][, area_szn := NULL]

message(' - aggiungo SZN senza OMI (vedere script `23b` dedicato) ')
ycna <- yz |> subset( !SZNid %in% rbindlist(list( yc1, yc2t[N==1] ))[, SZNid] )
ycna <- data.table(ycna |> st_drop_geometry(), yb[ycna |> st_nearest_feature(yb),] |> st_drop_geometry(), NA, NA, 0)

message(' - salvo file di ricerca SZN-OMI... ')
yc <- rbindlist(list( yc1, yc2t, ycna), use.names = FALSE) |> merge(masteRgeo::sezioni[, .(SZN, SZNid)], all.x = TRUE)
setcolorder(yc, 'SZN')
setorderv(yc, c('SZN', 'N'))
salva_dts_pkg(yc, 'szn_omi', dbn = 'omi')

message('\nDissolvo i confini Sezioni semplificati con associazioni...')
yb <-  st_read(file.path(bnd_path, '20.gpkg'), 'SZN', quiet = TRUE) |> subset(select = SZNid)
yc <- yc[N <= 1, .(SZNid, OMI)]
yb <- yb |> merge(yc) |> masteRgis::stm_dissolvi('OMI')

message('Salvo versione finale...')
salva_dts_pkg(yb, 'OMIz', as_db = FALSE, csv_in_pkg = FALSE)
st_write(yb, file.path(bnd_path, '20.gpkg'), 'OMIz', append = FALSE, quiet = TRUE)

message('\nElaborazioni Terminate. Pulizia Ambiente ed Uscita')
rm(list = ls())
gc()