Social Vulnerability Spain - Copernicus data

Contents

Social Vulnerability Spain - 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:"

'exactextractr'
'terra'
'sf'
'pacman'

### Output directory

# create the pipeline directory if it does not exist
pipeline_dir <- file.path("../..","2_pipeline","Spain","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 <- "Spain"

Process#

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

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

# census output areas
census_areas <- st_read("../../0_data/boundaries/Spain/census_sections/2021/Seccionado_2021/SECC_CE_20210101.shp")

Reading layer `SECC_CE_20210101' from data source 
  `/Cities/0_data/boundaries/Spain/census_sections/2021/Seccionado_2021/SECC_CE_20210101.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 36333 features and 18 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -1004502 ymin: 3132130 xmax: 1126932 ymax: 4859240
Projected CRS: ETRS89 / UTM zone 30N

# 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: IMD"
[1] "../../2_pipeline/Spain/1b_Copernicus/2021/Spain_IMD.tif already created"
[1] "Processing: TCD"
[1] "../../2_pipeline/Spain/1b_Copernicus/2021/Spain_TCD.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-03-13 18:36:33 GMT"

Warning message in .local(x, y, ...):
“Polygons transformed to raster CRS (EPSG:3035)”
Cannot preload entire working area of 34899458220 cells with max_cells_in_memory = 3e+07. Raster values will be read for each feature individually.

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

Warning message in .local(x, y, ...):
“Polygons transformed to raster CRS (EPSG:3035)”
Cannot preload entire working area of 34899458220 cells with max_cells_in_memory = 3e+07. Raster values will be read for each feature individually.

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

Registered S3 method overwritten by 'geojsonsf':
  method        from   
  print.geojson geojson

A sf: 6 × 23
	CUSEC	CUMUN	CSEC	CDIS	CMUN	CPRO	CCA	CUDIS	CLAU2	NPRO	geometry	⋯	CNUT2	CNUT3	ESTADO	OBS	NMUN	geometry	cell_area_m2	imp_area_m2	imp_percent	impervious
	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	⋯	<chr>	<chr>	<chr>	<chr>	<chr>	<MULTIPOLYGON [m]>	<dbl>	<dbl>	<dbl>	<dbl[,1]>
1	0100101001	01001	001	01	001	01	16	0100101	01001	Araba/Álava	MULTIPOLYGON (((539753 4743...	⋯	1	1	I	NA	Alegría-Dulantzi	MULTIPOLYGON (((539753 4743...	10819948.0	196633.8	1.817	-1.2163532
2	0100101002	01001	002	01	001	01	16	0100101	01001	Araba/Álava	MULTIPOLYGON (((539559.7 47...	⋯	1	1	I	NA	Alegría-Dulantzi	MULTIPOLYGON (((539559.7 47...	9120120.0	671134.9	7.359	-1.0566674
3	0100201001	01002	001	01	002	01	16	0100201	01002	Araba/Álava	MULTIPOLYGON (((503618.6 47...	⋯	1	1	I	NA	Amurrio	MULTIPOLYGON (((503618.6 47...	34817732.0	338359.4	0.972	-1.2407009
4	0100201002	01002	002	01	002	01	16	0100201	01002	Araba/Álava	MULTIPOLYGON (((505948.3 47...	⋯	1	1	I	NA	Amurrio	MULTIPOLYGON (((505948.3 47...	39335452.0	918302.3	2.335	-1.2014277
5	0100201003	01002	003	01	002	01	16	0100201	01002	Araba/Álava	MULTIPOLYGON (((499919.5 47...	⋯	1	1	I	NA	Amurrio	MULTIPOLYGON (((499919.5 47...	850154.1	97578.6	11.478	-0.9379836
6	0100201004	01002	004	01	002	01	16	0100201	01002	Araba/Álava	MULTIPOLYGON (((500127.9 47...	⋯	1	1	I	NA	Amurrio	MULTIPOLYGON (((500127.9 47...	492263.8	112533.7	22.860	-0.6100255

[1] "2025-03-14 02:01:09 GMT"

# 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/Spain/1b_Copernicus/2021/census_areas_IMD.geojson' using driver `GeoJSON'
Writing layer `census_areas_IMD' to data source 
  `../../2_pipeline/Spain/1b_Copernicus/2021/census_areas_IMD.geojson' using driver `GeoJSON'
Writing 36333 features with 22 fields and geometry type Multi 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-03-14 02:01:58 GMT"

Warning message in .local(x, y, ...):
“Polygons transformed to raster CRS (EPSG:3035)”
Cannot preload entire working area of 34899458220 cells with max_cells_in_memory = 3e+07. Raster values will be read for each feature individually.

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

Warning message in .local(x, y, ...):
“Polygons transformed to raster CRS (EPSG:3035)”
Cannot preload entire working area of 34899458220 cells with max_cells_in_memory = 3e+07. Raster values will be read for each feature individually.

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

A sf: 6 × 23
	CUSEC	CUMUN	CSEC	CDIS	CMUN	CPRO	CCA	CUDIS	CLAU2	NPRO	geometry	⋯	CNUT2	CNUT3	ESTADO	OBS	NMUN	geometry	cell_area_m2	tcd_area_m2	tcd_percent	tree_cover_density
	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	⋯	<chr>	<chr>	<chr>	<chr>	<chr>	<MULTIPOLYGON [m]>	<dbl>	<dbl>	<dbl>	<dbl[,1]>
1	0100101001	01001	001	01	001	01	16	0100101	01001	Araba/Álava	MULTIPOLYGON (((539753 4743...	⋯	1	1	I	NA	Alegría-Dulantzi	MULTIPOLYGON (((539753 4743...	10819948.0	3747061.00	34.631	-1.8991377
2	0100101002	01001	002	01	001	01	16	0100101	01001	Araba/Álava	MULTIPOLYGON (((539559.7 47...	⋯	1	1	I	NA	Alegría-Dulantzi	MULTIPOLYGON (((539559.7 47...	9120120.0	1123516.12	12.319	-0.1728744
3	0100201001	01002	001	01	002	01	16	0100201	01002	Araba/Álava	MULTIPOLYGON (((503618.6 47...	⋯	1	1	I	NA	Amurrio	MULTIPOLYGON (((503618.6 47...	34817732.0	12273660.00	35.251	-1.9471067
4	0100201002	01002	002	01	002	01	16	0100201	01002	Araba/Álava	MULTIPOLYGON (((505948.3 47...	⋯	1	1	I	NA	Amurrio	MULTIPOLYGON (((505948.3 47...	39335452.0	18191462.00	46.247	-2.7978593
5	0100201003	01002	003	01	002	01	16	0100201	01002	Araba/Álava	MULTIPOLYGON (((499919.5 47...	⋯	1	1	I	NA	Amurrio	MULTIPOLYGON (((499919.5 47...	850154.1	45916.07	5.401	0.3623663
6	0100201004	01002	004	01	002	01	16	0100201	01002	Araba/Álava	MULTIPOLYGON (((500127.9 47...	⋯	1	1	I	NA	Amurrio	MULTIPOLYGON (((500127.9 47...	492263.8	36810.55	7.478	0.2016703

[1] "2025-03-14 11:49:32 GMT"

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/Spain/1b_Copernicus/2021/census_areas_TCD.geojson' using driver `GeoJSON'
Writing layer `census_areas_TCD' to data source 
  `../../2_pipeline/Spain/1b_Copernicus/2021/census_areas_TCD.geojson' using driver `GeoJSON'
Writing 36333 features with 22 fields and geometry type Multi Polygon.

END