{ "cells": [ { "cell_type": "markdown", "id": "9065e87a-36f6-49f6-a3d2-49095355f564", "metadata": {}, "source": [ "# Social Vulnerability Italy - Copernicus data\n", "Processing of the Copernicus data into census output areas" ] }, { "cell_type": "markdown", "id": "17246fe7-bfca-4a51-b64d-2b337c9c239a", "metadata": {}, "source": [ "## Environment" ] }, { "cell_type": "markdown", "id": "a62367fa-e41d-47c0-be8b-b047093078e6", "metadata": {}, "source": [ "### R Libraries\n", "Any required R libraries are imported into the kernal:" ] }, { "cell_type": "code", "execution_count": 1, "id": "abb8d50e", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "Loading required package: pacman\n", "\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "[1] \"Loaded Packages:\"\n" ] }, { "data": { "text/html": [ "\n", "
  1. 'exactextractr'
  2. 'terra'
  3. 'sf'
  4. 'pacman'
\n" ], "text/latex": [ "\\begin{enumerate*}\n", "\\item 'exactextractr'\n", "\\item 'terra'\n", "\\item 'sf'\n", "\\item 'pacman'\n", "\\end{enumerate*}\n" ], "text/markdown": [ "1. 'exactextractr'\n", "2. 'terra'\n", "3. 'sf'\n", "4. 'pacman'\n", "\n", "\n" ], "text/plain": [ "[1] \"exactextractr\" \"terra\" \"sf\" \"pacman\" " ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Load R libraries\n", "if(!require(\"pacman\"))\n", " install.packages(\"pacman\")\n", "\n", "p_load(\"sf\", \"terra\", \"exactextractr\")\n", "\n", "print(\"Loaded Packages:\")\n", "p_loaded()" ] }, { "cell_type": "code", "execution_count": 2, "id": "00d29245-01af-4b94-9d8e-72e2a2c194a1", "metadata": {}, "outputs": [], "source": [ "### Output directory" ] }, { "cell_type": "code", "execution_count": 3, "id": "09702b65", "metadata": {}, "outputs": [], "source": [ "# create the pipeline directory if it does not exist\n", "pipeline_dir <- file.path(\"../..\",\"2_pipeline\",\"Italy\",\"Milan\",\"1b_Copernicus\",\"2021\")\n", "if(!dir.exists(pipeline_dir)){\n", " dir.create(pipeline_dir, recursive = TRUE)\n", " print(paste0(pipeline_dir, \" created\"))\n", "}\n", "\n", "# set country variable which is used within the output filename\n", "country <- \"Italy\"" ] }, { "cell_type": "markdown", "id": "7fb7048a-ed1b-4bbb-a43d-b103596ec774", "metadata": {}, "source": [ "## Process" ] }, { "cell_type": "code", "execution_count": 4, "id": "4cf607ac", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Reading layer `Sezioni_Censimento_2011' from data source \n", " `/Cities/0_data/boundaries/Italy/Milan/2021/ds98_infogeo_sezioni_censimento_localizzazione_2011c/Sezioni_Censimento_2011.shp' \n", " using driver `ESRI Shapefile'\n", "Simple feature collection with 6079 features and 11 fields\n", "Geometry type: POLYGON\n", "Dimension: XY\n", "Bounding box: xmin: 1503202 ymin: 5025952 xmax: 1521750 ymax: 5042528\n", "Projected CRS: Monte Mario / Italy zone 1\n" ] } ], "source": [ "# Copernicus data folder\n", "copernicus_data_path <- \"../../0_data/Copernicus/Italy\"\n", "\n", "# High resolution layers - IMPERVIOUSNESS, TREE COVER DENSITY\n", "datasets <- c(\"GRA\", \"IBU\",\"IMD\", \"TCD\", \"WAW\")\n", "\n", "# census output areas\n", "census_areas <- st_read(\"../../0_data/boundaries/Italy/Milan/2021/ds98_infogeo_sezioni_censimento_localizzazione_2011c/Sezioni_Censimento_2011.shp\")" ] }, { "cell_type": "code", "execution_count": 5, "id": "f0e29052", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[1] \"Processing: GRA\"\n", "[1] \"../../2_pipeline/Italy/Milan/1b_Copernicus/2021/Italy_GRA.tif already created\"\n", "[1] \"Processing: IBU\"\n", "[1] \"../../2_pipeline/Italy/Milan/1b_Copernicus/2021/Italy_IBU.tif already created\"\n", "[1] \"Processing: IMD\"\n", "[1] \"../../2_pipeline/Italy/Milan/1b_Copernicus/2021/Italy_IMD.tif already created\"\n", "[1] \"Processing: TCD\"\n", "[1] \"../../2_pipeline/Italy/Milan/1b_Copernicus/2021/Italy_TCD.tif already created\"\n", "[1] \"Processing: WAW\"\n", "[1] \"../../2_pipeline/Italy/Milan/1b_Copernicus/2021/Italy_WAW.tif already created\"\n" ] } ], "source": [ "# need to join copernicus tiles together\n", "for(dataset in datasets){\n", " print(paste0(\"Processing: \", dataset))\n", " flush.console()\n", "\n", " # export mosaic filepath\n", " mosaic_output_file <- paste0(pipeline_dir, \"/\", country, \"_\", dataset, \".tif\")\n", "\n", " # have the tiles been joined together already?\n", " is_mosaic_created <- file.exists(mosaic_output_file)\n", "\n", " if(is_mosaic_created){\n", " print(paste0(mosaic_output_file, \" already created\"))\n", " flush.console()\n", " } else {\n", " # Join the tiles together and export it\n", " # Get the files\n", " tiles <- list.files(copernicus_data_path, full.names = TRUE, recursive = TRUE, pattern=paste0(\".*\", dataset, \".*\\\\.tif$\")) \n", "\n", " # Create an image catalog\n", " ic <- terra::sprc(lapply(tiles, rast))\n", "\n", " # Mosic the tiles\n", " icMosaic <- mosaic(ic, filename = mosaic_output_file, fun=\"max\", overwrite=TRUE)\n", " #icMosaic <- merge(ic, filename = mosaic_output_file, overwrite=TRUE) # quicker\n", "\n", " print(paste0(mosaic_output_file, \" created\"))\n", " } \n", "}" ] }, { "cell_type": "code", "execution_count": 6, "id": "f7f1ac5b", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[1] \"2025-04-29 18:38:38 IST\"" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stderr", "output_type": "stream", "text": [ "Warning message in .local(x, y, ...):\n", "“Polygons transformed to raster CRS (EPSG:3035)”\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "\n", " |======================================================================| 100%" ] }, { "name": "stderr", "output_type": "stream", "text": [ "Warning message in .local(x, y, ...):\n", "“Polygons transformed to raster CRS (EPSG:3035)”\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "\n", " |======================================================================| 100%" ] }, { "name": "stderr", "output_type": "stream", "text": [ "Registered S3 method overwritten by 'geojsonsf':\n", " method from \n", " print.geojson geojson\n", "\n" ] }, { "data": { "application/geo+json": { "features": [ { "geometry": { "coordinates": [ [ [ 1514122.1494, 5034191.7775 ], [ 1514080.2244, 5034198.6974 ], [ 1514097.3244, 5034251.6125 ], [ 1514128.8944, 5034259.6974 ], [ 1514154.2894, 5034262.2775 ], [ 1514166.3394, 5034198.5825 ], [ 1514223.2944, 5034135.3425 ], [ 1514193.4394, 5034113.3375 ], [ 1514122.1494, 5034191.7775 ] ] ], "type": "Polygon" }, "properties": { "ACE": 1, "OBJECTID": 49010, "POP_2010": 70, "PRO_COM": 15146, "SEZ": 236, "SEZ2011": 151460000236, "SHAPE_AREA": 8061.4683, "SHAPE_LEN": 449.2267, "TIPO_LOC": 1, "cell_area_m2": 8067.7612, "imp_area_m2": 6153.9976, "imp_percent": 76.279, "impervious": [ [ 0.3807 ] ], "mappa2": 3, "mappa2bis": 5 }, "type": "Feature" }, { "geometry": { "coordinates": [ [ [ 1514166.3394, 5034198.5825 ], [ 1514154.2894, 5034262.2775 ], [ 1514216.5843, 5034267.4875 ], [ 1514217.1244, 5034212.3574 ], [ 1514223.2944, 5034135.3425 ], [ 1514166.3394, 5034198.5825 ] ] ], "type": "Polygon" }, "properties": { "ACE": 1, "OBJECTID": 49011, "POP_2010": 173, "PRO_COM": 15146, "SEZ": 237, "SEZ2011": 151460000237, "SHAPE_AREA": 5416.9115, "SHAPE_LEN": 344.8344, "TIPO_LOC": 1, "cell_area_m2": 5421.1401, "imp_area_m2": 4819.959, "imp_percent": 88.91, "impervious": [ [ 0.969 ] ], "mappa2": 12, "mappa2bis": 7.79 }, "type": "Feature" }, { "geometry": { "coordinates": [ [ [ 1514365.5094, 5034211.7375 ], [ 1514391.3794, 5034253.2325 ], [ 1514436.2193, 5034311.4374 ], [ 1514478.2044, 5034276.6075 ], [ 1514527.0094, 5034256.9874 ], [ 1514512.4793, 5034241.8324 ], [ 1514485.3093, 5034219.5725 ], [ 1514459.5993, 5034200.5175 ], [ 1514403.4244, 5034179.1675 ], [ 1514329.6044, 5034166.8524 ], [ 1514365.5094, 5034211.7375 ] ] ], "type": "Polygon" }, "properties": { "ACE": 1, "OBJECTID": 49015, "POP_2010": 160, "PRO_COM": 15146, "SEZ": 241, "SEZ2011": 151460000241, "SHAPE_AREA": 12107.8581, "SHAPE_LEN": 510.055, "TIPO_LOC": 1, "cell_area_m2": 12117.3086, "imp_area_m2": 11727.0703, "imp_percent": 96.779, "impervious": [ [ 1.3356 ] ], "mappa2": 13, "mappa2bis": 9.15 }, "type": "Feature" }, { "geometry": { "coordinates": [ [ [ 1514508.9844, 5034317.4075 ], [ 1514519.4943, 5034345.9125 ], [ 1514519.0194, 5034355.3075 ], [ 1514515.0144, 5034370.0975 ], [ 1514517.1093, 5034394.0074 ], [ 1514524.4343, 5034404.8324 ], [ 1514548.1443, 5034426.0825 ], [ 1514565.2844, 5034443.7924 ], [ 1514634.7144, 5034374.2824 ], [ 1514592.4194, 5034324.9674 ], [ 1514554.3694, 5034286.3525 ], [ 1514508.9844, 5034317.4075 ] ] ], "type": "Polygon" }, "properties": { "ACE": 1, "OBJECTID": 49018, "POP_2010": 105, "PRO_COM": 15146, "SEZ": 244, "SEZ2011": 151460000244, "SHAPE_AREA": 11178.7036, "SHAPE_LEN": 421.0807, "TIPO_LOC": 1, "cell_area_m2": 11187.4277, "imp_area_m2": 10128.7168, "imp_percent": 90.537, "impervious": [ [ 1.0448 ] ], "mappa2": 22, "mappa2bis": 27.16 }, "type": "Feature" }, { "geometry": { "coordinates": [ [ [ 1515504.8093, 5035162.9424 ], [ 1515522.4442, 5035177.1924 ], [ 1515587.0692, 5035092.3774 ], [ 1515562.8343, 5035076.2224 ], [ 1515504.8093, 5035162.9424 ] ] ], "type": "Polygon" }, "properties": { "ACE": 1, "OBJECTID": 49241, "POP_2010": 6, "PRO_COM": 15146, "SEZ": 151, "SEZ2011": 151460000151, "SHAPE_AREA": 2727.7693, "SHAPE_LEN": 262.7678, "TIPO_LOC": 1, "cell_area_m2": 2729.8962, "imp_area_m2": 1873.9333, "imp_percent": 68.645, "impervious": [ [ 0.0251 ] ], "mappa2": 1, "mappa2bis": 25 }, "type": "Feature" }, { "geometry": { "coordinates": [ [ [ 1515423.2093, 5034342.9425 ], [ 1515464.6642, 5034406.1175 ], [ 1515511.6093, 5034362.8924 ], [ 1515442.1493, 5034260.4574 ], [ 1515393.2542, 5034292.5975 ], [ 1515423.2093, 5034342.9425 ] ] ], "type": "Polygon" }, "properties": { "ACE": 1, "OBJECTID": 49253, "POP_2010": 0, "PRO_COM": 15146, "SEZ": 163, "SEZ2011": 151460000163, "SHAPE_AREA": 7925.5947, "SHAPE_LEN": 380.2311, "TIPO_LOC": 1, "cell_area_m2": 7931.7739, "imp_area_m2": 7690.0742, "imp_percent": 96.953, "impervious": [ [ 1.3437 ] ], "mappa2": 0, "mappa2bis": 0 }, "type": "Feature" } ], "type": "FeatureCollection" }, "text/html": [ "\n", "\n", "\n", "\t\n", "\t\n", "\n", "\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\n", "
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
\n" ], "text/latex": [ "A sf: 6 × 16\n", "\\begin{tabular}{r|llllllllllllllll}\n", " & OBJECTID & PRO\\_COM & SEZ & TIPO\\_LOC & SEZ2011 & SHAPE\\_AREA & SHAPE\\_LEN & POP\\_2010 & ACE & mappa2 & mappa2bis & geometry & cell\\_area\\_m2 & imp\\_area\\_m2 & imp\\_percent & impervious\\\\\n", " & & & & & & & & & & & & & & & & \\\\\n", "\\hline\n", "\t1 & 49010 & 15146 & 236 & 1 & 1.5146e+11 & 8061.468 & 449.2267 & 70 & 1 & 3 & 5.00 & POLYGON ((1514122 5034192, ... & 8067.761 & 6153.998 & 76.279 & 0.38065989\\\\\n", "\t2 & 49011 & 15146 & 237 & 1 & 1.5146e+11 & 5416.912 & 344.8344 & 173 & 1 & 12 & 7.79 & POLYGON ((1514166 5034199, ... & 5421.140 & 4819.959 & 88.910 & 0.96901330\\\\\n", "\t3 & 49015 & 15146 & 241 & 1 & 1.5146e+11 & 12107.858 & 510.0550 & 160 & 1 & 13 & 9.15 & POLYGON ((1514366 5034212, ... & 12117.309 & 11727.070 & 96.779 & 1.33555221\\\\\n", "\t4 & 49018 & 15146 & 244 & 1 & 1.5146e+11 & 11178.704 & 421.0807 & 105 & 1 & 22 & 27.16 & POLYGON ((1514509 5034317, ... & 11187.428 & 10128.717 & 90.537 & 1.04479914\\\\\n", "\t5 & 49241 & 15146 & 151 & 1 & 1.5146e+11 & 2727.769 & 262.7678 & 6 & 1 & 1 & 25.00 & POLYGON ((1515505 5035163, ... & 2729.896 & 1873.933 & 68.645 & 0.02506731\\\\\n", "\t6 & 49253 & 15146 & 163 & 1 & 1.5146e+11 & 7925.595 & 380.2311 & 0 & 1 & 0 & 0.00 & POLYGON ((1515423 5034343, ... & 7931.774 & 7690.074 & 96.953 & 1.34365715\\\\\n", "\\end{tabular}\n" ], "text/markdown": [ "\n", "A sf: 6 × 16\n", "\n", "| | OBJECTID <dbl> | PRO_COM <dbl> | SEZ <dbl> | TIPO_LOC <int> | SEZ2011 <dbl> | SHAPE_AREA <dbl> | SHAPE_LEN <dbl> | POP_2010 <dbl> | ACE <dbl> | mappa2 <dbl> | mappa2bis <dbl> | geometry <POLYGON [m]> | cell_area_m2 <dbl> | imp_area_m2 <dbl> | imp_percent <dbl> | impervious <dbl[,1]> |\n", "|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|\n", "| 1 | 49010 | 15146 | 236 | 1 | 1.5146e+11 | 8061.468 | 449.2267 | 70 | 1 | 3 | 5.00 | POLYGON ((1514122 5034192, ... | 8067.761 | 6153.998 | 76.279 | 0.38065989 |\n", "| 2 | 49011 | 15146 | 237 | 1 | 1.5146e+11 | 5416.912 | 344.8344 | 173 | 1 | 12 | 7.79 | POLYGON ((1514166 5034199, ... | 5421.140 | 4819.959 | 88.910 | 0.96901330 |\n", "| 3 | 49015 | 15146 | 241 | 1 | 1.5146e+11 | 12107.858 | 510.0550 | 160 | 1 | 13 | 9.15 | POLYGON ((1514366 5034212, ... | 12117.309 | 11727.070 | 96.779 | 1.33555221 |\n", "| 4 | 49018 | 15146 | 244 | 1 | 1.5146e+11 | 11178.704 | 421.0807 | 105 | 1 | 22 | 27.16 | POLYGON ((1514509 5034317, ... | 11187.428 | 10128.717 | 90.537 | 1.04479914 |\n", "| 5 | 49241 | 15146 | 151 | 1 | 1.5146e+11 | 2727.769 | 262.7678 | 6 | 1 | 1 | 25.00 | POLYGON ((1515505 5035163, ... | 2729.896 | 1873.933 | 68.645 | 0.02506731 |\n", "| 6 | 49253 | 15146 | 163 | 1 | 1.5146e+11 | 7925.595 | 380.2311 | 0 | 1 | 0 | 0.00 | POLYGON ((1515423 5034343, ... | 7931.774 | 7690.074 | 96.953 | 1.34365715 |\n", "\n" ], "text/plain": [ " OBJECTID PRO_COM SEZ TIPO_LOC SEZ2011 SHAPE_AREA SHAPE_LEN POP_2010 ACE\n", "1 49010 15146 236 1 1.5146e+11 8061.468 449.2267 70 1 \n", "2 49011 15146 237 1 1.5146e+11 5416.912 344.8344 173 1 \n", "3 49015 15146 241 1 1.5146e+11 12107.858 510.0550 160 1 \n", "4 49018 15146 244 1 1.5146e+11 11178.704 421.0807 105 1 \n", "5 49241 15146 151 1 1.5146e+11 2727.769 262.7678 6 1 \n", "6 49253 15146 163 1 1.5146e+11 7925.595 380.2311 0 1 \n", " mappa2 mappa2bis geometry cell_area_m2 imp_area_m2\n", "1 3 5.00 POLYGON ((1514122 5034192, ... 8067.761 6153.998 \n", "2 12 7.79 POLYGON ((1514166 5034199, ... 5421.140 4819.959 \n", "3 13 9.15 POLYGON ((1514366 5034212, ... 12117.309 11727.070 \n", "4 22 27.16 POLYGON ((1514509 5034317, ... 11187.428 10128.717 \n", "5 1 25.00 POLYGON ((1515505 5035163, ... 2729.896 1873.933 \n", "6 0 0.00 POLYGON ((1515423 5034343, ... 7931.774 7690.074 \n", " imp_percent impervious\n", "1 76.279 0.38065989\n", "2 88.910 0.96901330\n", "3 96.779 1.33555221\n", "4 90.537 1.04479914\n", "5 68.645 0.02506731\n", "6 96.953 1.34365715" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/plain": [ "[1] \"2025-04-29 18:39:09 IST\"" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "Sys.time()\n", "census_areas_impervious <- census_areas\n", "\n", "# Impervious surface\n", "copRaster <- rast(paste0(pipeline_dir, \"/\", country, \"_\", \"IMD\", \".tif\"))\n", "\n", "# calculate the area of the raster cell (m2)\n", "cellArea <- res(copRaster)[1]*res(copRaster)[2]\n", "\n", "# calculate the number of cells within each output area\n", "census_areas_impervious$cell_area_m2 <- exact_extract(copRaster, census_areas_impervious, 'count', coverage_area = TRUE, progress = TRUE)\n", "\n", "# calculate the area of each output area that is impervious\n", "census_areas_impervious$imp_area_m2 <- exact_extract(copRaster, census_areas_impervious, 'sum', progress = TRUE)\n", "\n", "# calculate the percentage of each output area that is impervious\n", "census_areas_impervious$imp_percent <- round((census_areas_impervious$imp_area_m2/census_areas_impervious$cell_area_m2)*100, 3)\n", "\n", "# calculate Z score\n", "census_areas_impervious$impervious <- scale(census_areas_impervious$imp_percent)\n", "\n", "head(census_areas_impervious)\n", "Sys.time()" ] }, { "cell_type": "code", "execution_count": 7, "id": "ecef30ee", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Deleting source `../../2_pipeline/Italy/Milan/1b_Copernicus/2021/census_areas_IMD.geojson' using driver `GeoJSON'\n", "Writing layer `census_areas_IMD' to data source \n", " `../../2_pipeline/Italy/Milan/1b_Copernicus/2021/census_areas_IMD.geojson' using driver `GeoJSON'\n", "Writing 6079 features with 15 fields and geometry type Polygon.\n" ] } ], "source": [ "# output the impervious data as a geojson\n", "st_write(census_areas_impervious, file.path(pipeline_dir, \"census_areas_IMD.geojson\"), delete_dsn = TRUE)" ] }, { "cell_type": "code", "execution_count": 8, "id": "bea93c90", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[1] \"2025-04-29 18:39:10 IST\"" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stderr", "output_type": "stream", "text": [ "Warning message in .local(x, y, ...):\n", "“Polygons transformed to raster CRS (EPSG:3035)”\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "\n", " |======================================================================| 100%" ] }, { "name": "stderr", "output_type": "stream", "text": [ "Warning message in .local(x, y, ...):\n", "“Polygons transformed to raster CRS (EPSG:3035)”\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "\n", " |======================================================================| 100%" ] }, { "data": { "application/geo+json": { "features": [ { "geometry": { "coordinates": [ [ [ 1514122.1494, 5034191.7775 ], [ 1514080.2244, 5034198.6974 ], [ 1514097.3244, 5034251.6125 ], [ 1514128.8944, 5034259.6974 ], [ 1514154.2894, 5034262.2775 ], [ 1514166.3394, 5034198.5825 ], [ 1514223.2944, 5034135.3425 ], [ 1514193.4394, 5034113.3375 ], [ 1514122.1494, 5034191.7775 ] ] ], "type": "Polygon" }, "properties": { "ACE": 1, "OBJECTID": 49010, "POP_2010": 70, "PRO_COM": 15146, "SEZ": 236, "SEZ2011": 151460000236, "SHAPE_AREA": 8061.4683, "SHAPE_LEN": 449.2267, "TIPO_LOC": 1, "cell_area_m2": 8067.7612, "mappa2": 3, "mappa2bis": 5, "tcd_area_m2": 0, "tcd_percent": 0, "tree_cover_density": [ [ 0.6144 ] ] }, "type": "Feature" }, { "geometry": { "coordinates": [ [ [ 1514166.3394, 5034198.5825 ], [ 1514154.2894, 5034262.2775 ], [ 1514216.5843, 5034267.4875 ], [ 1514217.1244, 5034212.3574 ], [ 1514223.2944, 5034135.3425 ], [ 1514166.3394, 5034198.5825 ] ] ], "type": "Polygon" }, "properties": { "ACE": 1, "OBJECTID": 49011, "POP_2010": 173, "PRO_COM": 15146, "SEZ": 237, "SEZ2011": 151460000237, "SHAPE_AREA": 5416.9115, "SHAPE_LEN": 344.8344, "TIPO_LOC": 1, "cell_area_m2": 5421.1401, "mappa2": 12, "mappa2bis": 7.79, "tcd_area_m2": 0, "tcd_percent": 0, "tree_cover_density": [ [ 0.6144 ] ] }, "type": "Feature" }, { "geometry": { "coordinates": [ [ [ 1514365.5094, 5034211.7375 ], [ 1514391.3794, 5034253.2325 ], [ 1514436.2193, 5034311.4374 ], [ 1514478.2044, 5034276.6075 ], [ 1514527.0094, 5034256.9874 ], [ 1514512.4793, 5034241.8324 ], [ 1514485.3093, 5034219.5725 ], [ 1514459.5993, 5034200.5175 ], [ 1514403.4244, 5034179.1675 ], [ 1514329.6044, 5034166.8524 ], [ 1514365.5094, 5034211.7375 ] ] ], "type": "Polygon" }, "properties": { "ACE": 1, "OBJECTID": 49015, "POP_2010": 160, "PRO_COM": 15146, "SEZ": 241, "SEZ2011": 151460000241, "SHAPE_AREA": 12107.8581, "SHAPE_LEN": 510.055, "TIPO_LOC": 1, "cell_area_m2": 12117.3086, "mappa2": 13, "mappa2bis": 9.15, "tcd_area_m2": 0, "tcd_percent": 0, "tree_cover_density": [ [ 0.6144 ] ] }, "type": "Feature" }, { "geometry": { "coordinates": [ [ [ 1514508.9844, 5034317.4075 ], [ 1514519.4943, 5034345.9125 ], [ 1514519.0194, 5034355.3075 ], [ 1514515.0144, 5034370.0975 ], [ 1514517.1093, 5034394.0074 ], [ 1514524.4343, 5034404.8324 ], [ 1514548.1443, 5034426.0825 ], [ 1514565.2844, 5034443.7924 ], [ 1514634.7144, 5034374.2824 ], [ 1514592.4194, 5034324.9674 ], [ 1514554.3694, 5034286.3525 ], [ 1514508.9844, 5034317.4075 ] ] ], "type": "Polygon" }, "properties": { "ACE": 1, "OBJECTID": 49018, "POP_2010": 105, "PRO_COM": 15146, "SEZ": 244, "SEZ2011": 151460000244, "SHAPE_AREA": 11178.7036, "SHAPE_LEN": 421.0807, "TIPO_LOC": 1, "cell_area_m2": 11187.4277, "mappa2": 22, "mappa2bis": 27.16, "tcd_area_m2": 0, "tcd_percent": 0, "tree_cover_density": [ [ 0.6144 ] ] }, "type": "Feature" }, { "geometry": { "coordinates": [ [ [ 1515504.8093, 5035162.9424 ], [ 1515522.4442, 5035177.1924 ], [ 1515587.0692, 5035092.3774 ], [ 1515562.8343, 5035076.2224 ], [ 1515504.8093, 5035162.9424 ] ] ], "type": "Polygon" }, "properties": { "ACE": 1, "OBJECTID": 49241, "POP_2010": 6, "PRO_COM": 15146, "SEZ": 151, "SEZ2011": 151460000151, "SHAPE_AREA": 2727.7693, "SHAPE_LEN": 262.7678, "TIPO_LOC": 1, "cell_area_m2": 2729.8962, "mappa2": 1, "mappa2bis": 25, "tcd_area_m2": 0, "tcd_percent": 0, "tree_cover_density": [ [ 0.6144 ] ] }, "type": "Feature" }, { "geometry": { "coordinates": [ [ [ 1515423.2093, 5034342.9425 ], [ 1515464.6642, 5034406.1175 ], [ 1515511.6093, 5034362.8924 ], [ 1515442.1493, 5034260.4574 ], [ 1515393.2542, 5034292.5975 ], [ 1515423.2093, 5034342.9425 ] ] ], "type": "Polygon" }, "properties": { "ACE": 1, "OBJECTID": 49253, "POP_2010": 0, "PRO_COM": 15146, "SEZ": 163, "SEZ2011": 151460000163, "SHAPE_AREA": 7925.5947, "SHAPE_LEN": 380.2311, "TIPO_LOC": 1, "cell_area_m2": 7931.7739, "mappa2": 0, "mappa2bis": 0, "tcd_area_m2": 0, "tcd_percent": 0, "tree_cover_density": [ [ 0.6144 ] ] }, "type": "Feature" } ], "type": "FeatureCollection" }, "text/html": [ "\n", "\n", "\n", "\t\n", "\t\n", "\n", "\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\n", "
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
\n" ], "text/latex": [ "A sf: 6 × 16\n", "\\begin{tabular}{r|llllllllllllllll}\n", " & OBJECTID & PRO\\_COM & SEZ & TIPO\\_LOC & SEZ2011 & SHAPE\\_AREA & SHAPE\\_LEN & POP\\_2010 & ACE & mappa2 & mappa2bis & geometry & cell\\_area\\_m2 & tcd\\_area\\_m2 & tcd\\_percent & tree\\_cover\\_density\\\\\n", " & & & & & & & & & & & & & & & & \\\\\n", "\\hline\n", "\t1 & 49010 & 15146 & 236 & 1 & 1.5146e+11 & 8061.468 & 449.2267 & 70 & 1 & 3 & 5.00 & POLYGON ((1514122 5034192, ... & 8067.761 & 0 & 0 & 0.6144226\\\\\n", "\t2 & 49011 & 15146 & 237 & 1 & 1.5146e+11 & 5416.912 & 344.8344 & 173 & 1 & 12 & 7.79 & POLYGON ((1514166 5034199, ... & 5421.140 & 0 & 0 & 0.6144226\\\\\n", "\t3 & 49015 & 15146 & 241 & 1 & 1.5146e+11 & 12107.858 & 510.0550 & 160 & 1 & 13 & 9.15 & POLYGON ((1514366 5034212, ... & 12117.309 & 0 & 0 & 0.6144226\\\\\n", "\t4 & 49018 & 15146 & 244 & 1 & 1.5146e+11 & 11178.704 & 421.0807 & 105 & 1 & 22 & 27.16 & POLYGON ((1514509 5034317, ... & 11187.428 & 0 & 0 & 0.6144226\\\\\n", "\t5 & 49241 & 15146 & 151 & 1 & 1.5146e+11 & 2727.769 & 262.7678 & 6 & 1 & 1 & 25.00 & POLYGON ((1515505 5035163, ... & 2729.896 & 0 & 0 & 0.6144226\\\\\n", "\t6 & 49253 & 15146 & 163 & 1 & 1.5146e+11 & 7925.595 & 380.2311 & 0 & 1 & 0 & 0.00 & POLYGON ((1515423 5034343, ... & 7931.774 & 0 & 0 & 0.6144226\\\\\n", "\\end{tabular}\n" ], "text/markdown": [ "\n", "A sf: 6 × 16\n", "\n", "| | OBJECTID <dbl> | PRO_COM <dbl> | SEZ <dbl> | TIPO_LOC <int> | SEZ2011 <dbl> | SHAPE_AREA <dbl> | SHAPE_LEN <dbl> | POP_2010 <dbl> | ACE <dbl> | mappa2 <dbl> | mappa2bis <dbl> | geometry <POLYGON [m]> | cell_area_m2 <dbl> | tcd_area_m2 <dbl> | tcd_percent <dbl> | tree_cover_density <dbl[,1]> |\n", "|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|\n", "| 1 | 49010 | 15146 | 236 | 1 | 1.5146e+11 | 8061.468 | 449.2267 | 70 | 1 | 3 | 5.00 | POLYGON ((1514122 5034192, ... | 8067.761 | 0 | 0 | 0.6144226 |\n", "| 2 | 49011 | 15146 | 237 | 1 | 1.5146e+11 | 5416.912 | 344.8344 | 173 | 1 | 12 | 7.79 | POLYGON ((1514166 5034199, ... | 5421.140 | 0 | 0 | 0.6144226 |\n", "| 3 | 49015 | 15146 | 241 | 1 | 1.5146e+11 | 12107.858 | 510.0550 | 160 | 1 | 13 | 9.15 | POLYGON ((1514366 5034212, ... | 12117.309 | 0 | 0 | 0.6144226 |\n", "| 4 | 49018 | 15146 | 244 | 1 | 1.5146e+11 | 11178.704 | 421.0807 | 105 | 1 | 22 | 27.16 | POLYGON ((1514509 5034317, ... | 11187.428 | 0 | 0 | 0.6144226 |\n", "| 5 | 49241 | 15146 | 151 | 1 | 1.5146e+11 | 2727.769 | 262.7678 | 6 | 1 | 1 | 25.00 | POLYGON ((1515505 5035163, ... | 2729.896 | 0 | 0 | 0.6144226 |\n", "| 6 | 49253 | 15146 | 163 | 1 | 1.5146e+11 | 7925.595 | 380.2311 | 0 | 1 | 0 | 0.00 | POLYGON ((1515423 5034343, ... | 7931.774 | 0 | 0 | 0.6144226 |\n", "\n" ], "text/plain": [ " OBJECTID PRO_COM SEZ TIPO_LOC SEZ2011 SHAPE_AREA SHAPE_LEN POP_2010 ACE\n", "1 49010 15146 236 1 1.5146e+11 8061.468 449.2267 70 1 \n", "2 49011 15146 237 1 1.5146e+11 5416.912 344.8344 173 1 \n", "3 49015 15146 241 1 1.5146e+11 12107.858 510.0550 160 1 \n", "4 49018 15146 244 1 1.5146e+11 11178.704 421.0807 105 1 \n", "5 49241 15146 151 1 1.5146e+11 2727.769 262.7678 6 1 \n", "6 49253 15146 163 1 1.5146e+11 7925.595 380.2311 0 1 \n", " mappa2 mappa2bis geometry cell_area_m2 tcd_area_m2\n", "1 3 5.00 POLYGON ((1514122 5034192, ... 8067.761 0 \n", "2 12 7.79 POLYGON ((1514166 5034199, ... 5421.140 0 \n", "3 13 9.15 POLYGON ((1514366 5034212, ... 12117.309 0 \n", "4 22 27.16 POLYGON ((1514509 5034317, ... 11187.428 0 \n", "5 1 25.00 POLYGON ((1515505 5035163, ... 2729.896 0 \n", "6 0 0.00 POLYGON ((1515423 5034343, ... 7931.774 0 \n", " tcd_percent tree_cover_density\n", "1 0 0.6144226 \n", "2 0 0.6144226 \n", "3 0 0.6144226 \n", "4 0 0.6144226 \n", "5 0 0.6144226 \n", "6 0 0.6144226 " ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/plain": [ "[1] \"2025-04-29 18:39:33 IST\"" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "Sys.time()\n", "census_areas_tree_cover <- census_areas\n", "\n", "# TREE COVER DENSITY\n", "# cell size = 100 m2, therefore a percentage value directly equals m2, i.e. 100% coverage = 100 m2 area of impervious surface\n", "copRaster <- rast(paste0(pipeline_dir, \"/\", country, \"_\", \"TCD\", \".tif\"))\n", "\n", "# calculate the area of the raster cell (m2)\n", "cellArea <- res(copRaster)[1]*res(copRaster)[2]\n", "\n", "# calculate the number of cells within each output area\n", "census_areas_tree_cover$cell_area_m2 <- exact_extract(copRaster, census_areas_tree_cover, 'count', coverage_area = TRUE, progress = TRUE)\n", "\n", "# calculate the area of each output area that has tree cover\n", "census_areas_tree_cover$tcd_area_m2 <- exact_extract(copRaster, census_areas_tree_cover, 'sum', progress = TRUE)\n", "\n", "# calculate the percentage of each output area that has tree cover\n", "census_areas_tree_cover$tcd_percent <- round((census_areas_tree_cover$tcd_area_m2/census_areas_tree_cover$cell_area_m2)*100, 3)\n", "\n", "# calculate Z score\n", "census_areas_tree_cover$tree_cover_density <- -scale(census_areas_tree_cover$tcd_percent)\n", "\n", "head(census_areas_tree_cover)\n", "Sys.time()" ] }, { "cell_type": "markdown", "id": "4c7b45ca-98c0-4fd7-ace7-db0c526ca01e", "metadata": {}, "source": [ "## Export" ] }, { "cell_type": "code", "execution_count": 9, "id": "f0bc25eb", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Deleting source `../../2_pipeline/Italy/Milan/1b_Copernicus/2021/census_areas_TCD.geojson' using driver `GeoJSON'\n", "Writing layer `census_areas_TCD' to data source \n", " `../../2_pipeline/Italy/Milan/1b_Copernicus/2021/census_areas_TCD.geojson' using driver `GeoJSON'\n", "Writing 6079 features with 15 fields and geometry type Polygon.\n" ] } ], "source": [ "# output the tree cover density data as a geojson\n", "st_write(census_areas_tree_cover, file.path(pipeline_dir, \"census_areas_TCD.geojson\"), delete_dsn = TRUE)" ] }, { "cell_type": "markdown", "id": "90b4c24a-8695-4340-b3bc-46b63dabc967", "metadata": {}, "source": [ "**END**" ] } ], "metadata": { "kernelspec": { "display_name": "R", "language": "R", "name": "ir" }, "language_info": { "codemirror_mode": "r", "file_extension": ".r", "mimetype": "text/x-r-source", "name": "R", "pygments_lexer": "r", "version": "4.3.3" } }, "nbformat": 4, "nbformat_minor": 5 }