Social Vulnerability Italy - Copernicus data

Social Vulnerability Italy - Copernicus data#

Processing of the Copernicus data into census output areas

Environment#

R Libraries#

Any required R libraries are imported into the kernal:

# Load R libraries
if(!require("pacman"))
    install.packages("pacman")

p_load("sf", "terra", "exactextractr")

print("Loaded Packages:")
p_loaded()

Loading required package: pacman

[1] "Loaded Packages:"
  1. 'exactextractr'
  2. 'terra'
  3. 'sf'
  4. 'pacman'
### Output directory

# create the pipeline directory if it does not exist
pipeline_dir <- file.path("../..","2_pipeline","Italy","Milan","1b_Copernicus","2021")
if(!dir.exists(pipeline_dir)){
    dir.create(pipeline_dir, recursive = TRUE)
    print(paste0(pipeline_dir, " created"))
}

# set country variable which is used within the output filename
country <- "Italy"

Process#

# Copernicus data folder
copernicus_data_path <- "../../0_data/Copernicus/Italy"

# High resolution layers - IMPERVIOUSNESS, TREE COVER DENSITY
datasets <- c("GRA", "IBU","IMD", "TCD", "WAW")

# census output areas
census_areas <- st_read("../../0_data/boundaries/Italy/Milan/2021/ds98_infogeo_sezioni_censimento_localizzazione_2011c/Sezioni_Censimento_2011.shp")

Reading layer `Sezioni_Censimento_2011' from data source 
  `/Cities/0_data/boundaries/Italy/Milan/2021/ds98_infogeo_sezioni_censimento_localizzazione_2011c/Sezioni_Censimento_2011.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 6079 features and 11 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 1503202 ymin: 5025952 xmax: 1521750 ymax: 5042528
Projected CRS: Monte Mario / Italy zone 1

# need to join copernicus tiles together
for(dataset in datasets){
    print(paste0("Processing: ", dataset))
    flush.console()

    # export mosaic filepath
    mosaic_output_file <- paste0(pipeline_dir, "/", country, "_", dataset, ".tif")

    # have the tiles been joined together already?
    is_mosaic_created <- file.exists(mosaic_output_file)

    if(is_mosaic_created){
        print(paste0(mosaic_output_file, " already created"))
        flush.console()
    } else {
        # Join the tiles together and export it
        # Get the files
        tiles <- list.files(copernicus_data_path, full.names = TRUE, recursive = TRUE, pattern=paste0(".*", dataset, ".*\\.tif$"))  

        # Create an image catalog
        ic <- terra::sprc(lapply(tiles, rast))

        # Mosic the tiles
        icMosaic <- mosaic(ic, filename = mosaic_output_file, fun="max", overwrite=TRUE)
        #icMosaic <- merge(ic, filename = mosaic_output_file, overwrite=TRUE) # quicker

        print(paste0(mosaic_output_file, " created"))
    }   
}

[1] "Processing: GRA"
[1] "../../2_pipeline/Italy/Milan/1b_Copernicus/2021/Italy_GRA.tif already created"
[1] "Processing: IBU"
[1] "../../2_pipeline/Italy/Milan/1b_Copernicus/2021/Italy_IBU.tif already created"
[1] "Processing: IMD"
[1] "../../2_pipeline/Italy/Milan/1b_Copernicus/2021/Italy_IMD.tif already created"
[1] "Processing: TCD"
[1] "../../2_pipeline/Italy/Milan/1b_Copernicus/2021/Italy_TCD.tif already created"
[1] "Processing: WAW"
[1] "../../2_pipeline/Italy/Milan/1b_Copernicus/2021/Italy_WAW.tif already created"

Sys.time()
census_areas_impervious <- census_areas

# Impervious surface
copRaster <- rast(paste0(pipeline_dir, "/", country, "_", "IMD", ".tif"))

# calculate the area of the raster cell (m2)
cellArea <- res(copRaster)[1]*res(copRaster)[2]

# calculate the number of cells within each output area
census_areas_impervious$cell_area_m2 <- exact_extract(copRaster, census_areas_impervious, 'count', coverage_area = TRUE, progress = TRUE)

# calculate the area of each output area that is impervious
census_areas_impervious$imp_area_m2 <- exact_extract(copRaster, census_areas_impervious, 'sum', progress = TRUE)

# calculate the percentage of each output area that is impervious
census_areas_impervious$imp_percent <- round((census_areas_impervious$imp_area_m2/census_areas_impervious$cell_area_m2)*100, 3)

# calculate Z score
census_areas_impervious$impervious <- scale(census_areas_impervious$imp_percent)

head(census_areas_impervious)
Sys.time()

[1] "2025-04-29 18:38:38 IST"

Warning message in .local(x, y, ...):
“Polygons transformed to raster CRS (EPSG:3035)”

 |======================================================================| 100%

Warning message in .local(x, y, ...):
“Polygons transformed to raster CRS (EPSG:3035)”

 |======================================================================| 100%

Registered S3 method overwritten by 'geojsonsf':
  method        from   
  print.geojson geojson
A sf: 6 × 16
OBJECTIDPRO_COMSEZTIPO_LOCSEZ2011SHAPE_AREASHAPE_LENPOP_2010ACEmappa2mappa2bisgeometrycell_area_m2imp_area_m2imp_percentimpervious
<dbl><dbl><dbl><int><dbl><dbl><dbl><dbl><dbl><dbl><dbl><POLYGON [m]><dbl><dbl><dbl><dbl[,1]>
1490101514623611.5146e+11 8061.468449.2267 701 3 5.00POLYGON ((1514122 5034192, ... 8067.761 6153.99876.2790.38065989
2490111514623711.5146e+11 5416.912344.8344173112 7.79POLYGON ((1514166 5034199, ... 5421.140 4819.95988.9100.96901330
3490151514624111.5146e+1112107.858510.0550160113 9.15POLYGON ((1514366 5034212, ...12117.30911727.07096.7791.33555221
4490181514624411.5146e+1111178.704421.080710512227.16POLYGON ((1514509 5034317, ...11187.42810128.71790.5371.04479914
5492411514615111.5146e+11 2727.769262.7678 61 125.00POLYGON ((1515505 5035163, ... 2729.896 1873.93368.6450.02506731
6492531514616311.5146e+11 7925.595380.2311 01 0 0.00POLYGON ((1515423 5034343, ... 7931.774 7690.07496.9531.34365715
 
[1] "2025-04-29 18:39:09 IST"

# output the impervious data as a geojson
st_write(census_areas_impervious, file.path(pipeline_dir, "census_areas_IMD.geojson"), delete_dsn = TRUE)

Deleting source `../../2_pipeline/Italy/Milan/1b_Copernicus/2021/census_areas_IMD.geojson' using driver `GeoJSON'
Writing layer `census_areas_IMD' to data source 
  `../../2_pipeline/Italy/Milan/1b_Copernicus/2021/census_areas_IMD.geojson' using driver `GeoJSON'
Writing 6079 features with 15 fields and geometry type Polygon.

Sys.time()
census_areas_tree_cover <- census_areas

# TREE COVER DENSITY
# cell size = 100 m2, therefore a percentage value directly equals m2, i.e. 100% coverage = 100 m2 area of impervious surface
copRaster <- rast(paste0(pipeline_dir, "/", country, "_", "TCD", ".tif"))

# calculate the area of the raster cell (m2)
cellArea <- res(copRaster)[1]*res(copRaster)[2]

# calculate the number of cells within each output area
census_areas_tree_cover$cell_area_m2 <- exact_extract(copRaster, census_areas_tree_cover, 'count', coverage_area = TRUE, progress = TRUE)

# calculate the area of each output area that has tree cover
census_areas_tree_cover$tcd_area_m2 <- exact_extract(copRaster, census_areas_tree_cover, 'sum', progress = TRUE)

# calculate the percentage of each output area that has tree cover
census_areas_tree_cover$tcd_percent <- round((census_areas_tree_cover$tcd_area_m2/census_areas_tree_cover$cell_area_m2)*100, 3)

# calculate Z score
census_areas_tree_cover$tree_cover_density <- -scale(census_areas_tree_cover$tcd_percent)

head(census_areas_tree_cover)
Sys.time()

[1] "2025-04-29 18:39:10 IST"

Warning message in .local(x, y, ...):
“Polygons transformed to raster CRS (EPSG:3035)”

 |======================================================================| 100%

Warning message in .local(x, y, ...):
“Polygons transformed to raster CRS (EPSG:3035)”

 |======================================================================| 100%
A sf: 6 × 16
OBJECTIDPRO_COMSEZTIPO_LOCSEZ2011SHAPE_AREASHAPE_LENPOP_2010ACEmappa2mappa2bisgeometrycell_area_m2tcd_area_m2tcd_percenttree_cover_density
<dbl><dbl><dbl><int><dbl><dbl><dbl><dbl><dbl><dbl><dbl><POLYGON [m]><dbl><dbl><dbl><dbl[,1]>
1490101514623611.5146e+11 8061.468449.2267 701 3 5.00POLYGON ((1514122 5034192, ... 8067.761000.6144226
2490111514623711.5146e+11 5416.912344.8344173112 7.79POLYGON ((1514166 5034199, ... 5421.140000.6144226
3490151514624111.5146e+1112107.858510.0550160113 9.15POLYGON ((1514366 5034212, ...12117.309000.6144226
4490181514624411.5146e+1111178.704421.080710512227.16POLYGON ((1514509 5034317, ...11187.428000.6144226
5492411514615111.5146e+11 2727.769262.7678 61 125.00POLYGON ((1515505 5035163, ... 2729.896000.6144226
6492531514616311.5146e+11 7925.595380.2311 01 0 0.00POLYGON ((1515423 5034343, ... 7931.774000.6144226
 
[1] "2025-04-29 18:39:33 IST"

Export#

# output the tree cover density data as a geojson
st_write(census_areas_tree_cover, file.path(pipeline_dir, "census_areas_TCD.geojson"), delete_dsn = TRUE)

Deleting source `../../2_pipeline/Italy/Milan/1b_Copernicus/2021/census_areas_TCD.geojson' using driver `GeoJSON'
Writing layer `census_areas_TCD' to data source 
  `../../2_pipeline/Italy/Milan/1b_Copernicus/2021/census_areas_TCD.geojson' using driver `GeoJSON'
Writing 6079 features with 15 fields and geometry type Polygon.

END

Contents