Spaces:
Running
Running
Commit
·
607b7f5
1
Parent(s):
432797d
rounding
Browse files- get_zonal_stats.ipynb +0 -1056
get_zonal_stats.ipynb
DELETED
|
@@ -1,1056 +0,0 @@
|
|
| 1 |
-
{
|
| 2 |
-
"cells": [
|
| 3 |
-
{
|
| 4 |
-
"cell_type": "markdown",
|
| 5 |
-
"id": "39bf1de3-cba6-475a-a988-ad48e5af4a04",
|
| 6 |
-
"metadata": {},
|
| 7 |
-
"source": [
|
| 8 |
-
"# Get zonal stats "
|
| 9 |
-
]
|
| 10 |
-
},
|
| 11 |
-
{
|
| 12 |
-
"cell_type": "code",
|
| 13 |
-
"execution_count": null,
|
| 14 |
-
"id": "ba047a55-642d-4c27-a367-5f35f4406218",
|
| 15 |
-
"metadata": {},
|
| 16 |
-
"outputs": [],
|
| 17 |
-
"source": [
|
| 18 |
-
"import ibis\n",
|
| 19 |
-
"import ibis.selectors as s\n",
|
| 20 |
-
"from ibis import _\n",
|
| 21 |
-
"import fiona\n",
|
| 22 |
-
"import geopandas as gpd\n",
|
| 23 |
-
"import rioxarray\n",
|
| 24 |
-
"from shapely.geometry import box\n",
|
| 25 |
-
"\n",
|
| 26 |
-
"import rasterio\n",
|
| 27 |
-
"from rasterio.mask import mask\n",
|
| 28 |
-
"from rasterstats import zonal_stats\n",
|
| 29 |
-
"import pandas as pd\n",
|
| 30 |
-
"from joblib import Parallel, delayed\n",
|
| 31 |
-
"\n",
|
| 32 |
-
"con = ibis.duckdb.connect()\n",
|
| 33 |
-
"con.load_extension(\"spatial\")\n",
|
| 34 |
-
"threads = -1"
|
| 35 |
-
]
|
| 36 |
-
},
|
| 37 |
-
{
|
| 38 |
-
"cell_type": "code",
|
| 39 |
-
"execution_count": null,
|
| 40 |
-
"id": "8b5656db-2d1d-4ca8-826d-7588126e52e8",
|
| 41 |
-
"metadata": {},
|
| 42 |
-
"outputs": [],
|
| 43 |
-
"source": [
|
| 44 |
-
"# cropping US data to only CA \n",
|
| 45 |
-
"def crop_raster_to_bounds(tif_file, vector_gdf):\n",
|
| 46 |
-
" with rasterio.open(tif_file) as src:\n",
|
| 47 |
-
" # Get California's bounding box in the same CRS as the raster\n",
|
| 48 |
-
" california_bounds = vector_gdf.total_bounds\n",
|
| 49 |
-
" california_bounds = rasterio.coords.BoundingBox(\n",
|
| 50 |
-
" *california_bounds\n",
|
| 51 |
-
" )\n",
|
| 52 |
-
" # Crop the raster to the California bounding box\n",
|
| 53 |
-
" out_image, out_transform = mask(src, [california_bounds], crop=True)\n",
|
| 54 |
-
" out_meta = src.meta.copy()\n",
|
| 55 |
-
" out_meta.update({\n",
|
| 56 |
-
" \"driver\": \"GTiff\",\n",
|
| 57 |
-
" \"height\": out_image.shape[1],\n",
|
| 58 |
-
" \"width\": out_image.shape[2],\n",
|
| 59 |
-
" \"transform\": out_transform\n",
|
| 60 |
-
" })\n",
|
| 61 |
-
" print(\"Unique values in cropped raster:\", np.unique(out_image))\n",
|
| 62 |
-
"\n",
|
| 63 |
-
" return out_image, out_meta\n"
|
| 64 |
-
]
|
| 65 |
-
},
|
| 66 |
-
{
|
| 67 |
-
"cell_type": "code",
|
| 68 |
-
"execution_count": null,
|
| 69 |
-
"id": "9a0e3446-16ac-40b0-9e34-db0157038c5a",
|
| 70 |
-
"metadata": {},
|
| 71 |
-
"outputs": [],
|
| 72 |
-
"source": [
|
| 73 |
-
"def big_zonal_stats(vec_file, tif_file, stats, col_name, n_jobs, verbose=10, timeout=10000):\n",
|
| 74 |
-
" gdf = gpd.read_parquet(vec_file)\n",
|
| 75 |
-
" if gdf.crs is None:\n",
|
| 76 |
-
" gdf = gdf.set_crs(\"EPSG:4326\")\n",
|
| 77 |
-
" gdf = gdf.rename(columns={\"geom\": \"geometry\"})\n",
|
| 78 |
-
" gdf = gdf.set_geometry(\"geometry\")\n",
|
| 79 |
-
" gdf = gdf[gdf[\"geometry\"].notna()].copy()\n",
|
| 80 |
-
"\n",
|
| 81 |
-
" with rasterio.open(tif_file) as src:\n",
|
| 82 |
-
" raster_crs = src.crs\n",
|
| 83 |
-
" gdf = gdf.to_crs(raster_crs) # Transform vector to raster CRS\n",
|
| 84 |
-
" \n",
|
| 85 |
-
" # CA bounding box + convert it to a polygon in raster CRS\n",
|
| 86 |
-
" california_polygon = box(*gdf.total_bounds)\n",
|
| 87 |
-
" \n",
|
| 88 |
-
" out_image, out_transform = mask(src, [california_polygon], crop=True, nodata=src.nodata)\n",
|
| 89 |
-
"\n",
|
| 90 |
-
" # If raster is 3D, select the first band\n",
|
| 91 |
-
" if out_image.ndim == 3:\n",
|
| 92 |
-
" out_image = out_image[0]\n",
|
| 93 |
-
"\n",
|
| 94 |
-
" # compute zonal statistics for each geometry slice\n",
|
| 95 |
-
" def get_stats(geom_slice):\n",
|
| 96 |
-
" geom = [geom_slice.geometry]\n",
|
| 97 |
-
" stats_result = zonal_stats(\n",
|
| 98 |
-
" geom, out_image, stats=stats, affine=out_transform, all_touched=True, nodata=src.nodata\n",
|
| 99 |
-
" )\n",
|
| 100 |
-
" return stats_result[0] if stats_result and stats_result[0].get(\"mean\") is not None else {'mean': None}\n",
|
| 101 |
-
"\n",
|
| 102 |
-
" output = [get_stats(row) for row in gdf.itertuples()]\n",
|
| 103 |
-
" gdf[col_name] = [res['mean'] for res in output]\n",
|
| 104 |
-
"\n",
|
| 105 |
-
" return gdf"
|
| 106 |
-
]
|
| 107 |
-
},
|
| 108 |
-
{
|
| 109 |
-
"cell_type": "code",
|
| 110 |
-
"execution_count": null,
|
| 111 |
-
"id": "ce66bae6-bac5-4837-9b01-fde16a00c303",
|
| 112 |
-
"metadata": {},
|
| 113 |
-
"outputs": [],
|
| 114 |
-
"source": [
|
| 115 |
-
"# getting local copies of data \n",
|
| 116 |
-
"# aws s3 cp s3://vizzuality/hfp-100/hfp_2021_100m_v1-2_cog.tif . --endpoint-url=https://data.source.coop\n",
|
| 117 |
-
"# aws s3 cp s3://vizzuality/lg-land-carbon-data/natcrop_bii_100m_cog.tif . --endpoint-url=https://data.source.coop\n",
|
| 118 |
-
"# aws s3 cp s3://vizzuality/lg-land-carbon-data/natcrop_fii_100m_cog.tif . --endpoint-url=https://data.source.coop\n",
|
| 119 |
-
"# aws s3 cp s3://vizzuality/lg-land-carbon-data/natcrop_expansion_100m_cog.tif . --endpoint-url=https://data.source.coop\n",
|
| 120 |
-
"# aws s3 cp s3://vizzuality/lg-land-carbon-data/natcrop_reduction_100m_cog.tif . --endpoint-url=https://data.source.coop\n",
|
| 121 |
-
"# aws s3 cp s3://cboettig/carbon/cogs/irrecoverable_c_total_2018.tif . --endpoint-url=https://data.source.coop\n",
|
| 122 |
-
"# aws s3 cp s3://cboettig/carbon/cogs/manageable_c_total_2018.tif . --endpoint-url=https://data.source.coop\n",
|
| 123 |
-
"# ! aws s3 cp s3://cboettig/justice40/disadvantaged-communities.parquet . --endpoint-url=https://data.source.coop\n",
|
| 124 |
-
"# minio/shared-biodiversity/redlist/cog/combined_sr_2022.tif\n",
|
| 125 |
-
"# /home/rstudio/minio/shared-biodiversity/redlist/cog/combined_rwr_2022.tif\n",
|
| 126 |
-
"# ! aws s3 cp s3://cboettig/social-vulnerability/svi2020_us_tract.parquet . --endpoint-url=https://data.source.coop\n"
|
| 127 |
-
]
|
| 128 |
-
},
|
| 129 |
-
{
|
| 130 |
-
"cell_type": "markdown",
|
| 131 |
-
"id": "531e7f88-1ce1-4027-b0ab-aab597e9a2b2",
|
| 132 |
-
"metadata": {},
|
| 133 |
-
"source": [
|
| 134 |
-
"# Biodiversity Data"
|
| 135 |
-
]
|
| 136 |
-
},
|
| 137 |
-
{
|
| 138 |
-
"cell_type": "code",
|
| 139 |
-
"execution_count": null,
|
| 140 |
-
"id": "66dec912-ad8a-41cf-a5c2-6ec9cc350984",
|
| 141 |
-
"metadata": {},
|
| 142 |
-
"outputs": [],
|
| 143 |
-
"source": [
|
| 144 |
-
"%%time\n",
|
| 145 |
-
"tif_file = 'SpeciesRichness_All.tif'\n",
|
| 146 |
-
"vec_file = \"/home/rstudio/github/ca-30x30/ca2024-30m.parquet\"\n",
|
| 147 |
-
"df = big_zonal_stats(vec_file, tif_file, stats = ['mean'], col_name = \"richness\", n_jobs=threads, verbose=0).to_parquet(\"cpad-stats-temp.parquet\")\n"
|
| 148 |
-
]
|
| 149 |
-
},
|
| 150 |
-
{
|
| 151 |
-
"cell_type": "code",
|
| 152 |
-
"execution_count": null,
|
| 153 |
-
"id": "b081ec1a-ea91-485e-95f9-12cd06c2002a",
|
| 154 |
-
"metadata": {},
|
| 155 |
-
"outputs": [],
|
| 156 |
-
"source": [
|
| 157 |
-
"%%time\n",
|
| 158 |
-
"tif_file = 'RSR_All.tif'\n",
|
| 159 |
-
"vec_file = './cpad-stats-temp.parquet'\n",
|
| 160 |
-
"df = big_zonal_stats(vec_file, tif_file, stats = ['mean'],\n",
|
| 161 |
-
" col_name = \"rsr\", n_jobs=threads, verbose=0).to_parquet(\"cpad-stats-temp.parquet\")"
|
| 162 |
-
]
|
| 163 |
-
},
|
| 164 |
-
{
|
| 165 |
-
"cell_type": "code",
|
| 166 |
-
"execution_count": null,
|
| 167 |
-
"id": "d5133f36-404e-4f6a-a90b-eb5f098e6f06",
|
| 168 |
-
"metadata": {},
|
| 169 |
-
"outputs": [],
|
| 170 |
-
"source": [
|
| 171 |
-
"%%time\n",
|
| 172 |
-
"tif_file = 'combined_sr_2022.tif'\n",
|
| 173 |
-
"vec_file = './cpad-stats-temp.parquet'\n",
|
| 174 |
-
"df = big_zonal_stats(vec_file, tif_file, stats = ['mean'], col_name = \"all_species_richness\", n_jobs=threads, verbose=0).to_parquet(\"cpad-stats-temp.parquet\")\n"
|
| 175 |
-
]
|
| 176 |
-
},
|
| 177 |
-
{
|
| 178 |
-
"cell_type": "code",
|
| 179 |
-
"execution_count": null,
|
| 180 |
-
"id": "2ce56a66-34e3-4f61-95ae-65d1f06bc468",
|
| 181 |
-
"metadata": {},
|
| 182 |
-
"outputs": [],
|
| 183 |
-
"source": [
|
| 184 |
-
"%%time\n",
|
| 185 |
-
"tif_file = 'combined_rwr_2022.tif'\n",
|
| 186 |
-
"vec_file = './cpad-stats-temp.parquet'\n",
|
| 187 |
-
"df = big_zonal_stats(vec_file, tif_file, stats = ['mean'], col_name = \"all_species_rwr\", n_jobs=threads, verbose=0).to_parquet(\"cpad-stats-temp.parquet\")\n"
|
| 188 |
-
]
|
| 189 |
-
},
|
| 190 |
-
{
|
| 191 |
-
"cell_type": "markdown",
|
| 192 |
-
"id": "6c129894-3775-4842-8767-f81a8f626d2c",
|
| 193 |
-
"metadata": {},
|
| 194 |
-
"source": [
|
| 195 |
-
"# Carbon Data"
|
| 196 |
-
]
|
| 197 |
-
},
|
| 198 |
-
{
|
| 199 |
-
"cell_type": "code",
|
| 200 |
-
"execution_count": null,
|
| 201 |
-
"id": "19c3e402-8712-450f-b3dd-af9d0c01689c",
|
| 202 |
-
"metadata": {},
|
| 203 |
-
"outputs": [],
|
| 204 |
-
"source": [
|
| 205 |
-
"%%time\n",
|
| 206 |
-
"tif_file = 'irrecoverable_c_total_2018.tif'\n",
|
| 207 |
-
"vec_file = './cpad-stats-temp.parquet'\n",
|
| 208 |
-
"df = big_zonal_stats(vec_file, tif_file, stats = ['mean'], col_name = \"irrecoverable_carbon\", n_jobs=threads, verbose=0).to_parquet(\"cpad-stats-temp.parquet\")\n",
|
| 209 |
-
"\n"
|
| 210 |
-
]
|
| 211 |
-
},
|
| 212 |
-
{
|
| 213 |
-
"cell_type": "code",
|
| 214 |
-
"execution_count": null,
|
| 215 |
-
"id": "c55c777a-48ce-4403-a171-cfc0d2351df6",
|
| 216 |
-
"metadata": {},
|
| 217 |
-
"outputs": [],
|
| 218 |
-
"source": [
|
| 219 |
-
"%%time\n",
|
| 220 |
-
"tif_file = 'manageable_c_total_2018.tif'\n",
|
| 221 |
-
"vec_file = './cpad-stats-temp.parquet'\n",
|
| 222 |
-
"df = big_zonal_stats(vec_file, tif_file, stats = ['mean'], col_name = \"manageable_carbon\", n_jobs=threads, verbose=0).to_parquet(\"cpad-stats-temp.parquet\")\n"
|
| 223 |
-
]
|
| 224 |
-
},
|
| 225 |
-
{
|
| 226 |
-
"cell_type": "code",
|
| 227 |
-
"execution_count": null,
|
| 228 |
-
"id": "33ac0fb7-2cde-448d-a634-1973e34ac14f",
|
| 229 |
-
"metadata": {},
|
| 230 |
-
"outputs": [],
|
| 231 |
-
"source": [
|
| 232 |
-
"%%time\n",
|
| 233 |
-
"tif_file = 'deforest_carbon_100m_cog.tif'\n",
|
| 234 |
-
"vec_file = './cpad-stats-temp.parquet'\n",
|
| 235 |
-
"df = big_zonal_stats(vec_file, tif_file, stats = ['mean'], \n",
|
| 236 |
-
" col_name = \"deforest_carbon\", n_jobs=threads, verbose=0).to_parquet(\"cpad-stats-temp.parquet\")\n"
|
| 237 |
-
]
|
| 238 |
-
},
|
| 239 |
-
{
|
| 240 |
-
"cell_type": "markdown",
|
| 241 |
-
"id": "096c00a8-57af-41d7-93cc-85d85414aa4f",
|
| 242 |
-
"metadata": {},
|
| 243 |
-
"source": [
|
| 244 |
-
"# Human Impact Data"
|
| 245 |
-
]
|
| 246 |
-
},
|
| 247 |
-
{
|
| 248 |
-
"cell_type": "code",
|
| 249 |
-
"execution_count": null,
|
| 250 |
-
"id": "d2a8c10f-e94b-4eef-940f-2af9599edee1",
|
| 251 |
-
"metadata": {},
|
| 252 |
-
"outputs": [],
|
| 253 |
-
"source": [
|
| 254 |
-
"%%time\n",
|
| 255 |
-
"tif_file = 'natcrop_bii_100m_cog.tif'\n",
|
| 256 |
-
"vec_file = './cpad-stats-temp.parquet'\n",
|
| 257 |
-
"df = big_zonal_stats(vec_file, tif_file, stats = ['mean'], \n",
|
| 258 |
-
" col_name = \"biodiversity_intactness_loss\", n_jobs=threads, verbose=0).to_parquet(\"cpad-stats-temp.parquet\")\n"
|
| 259 |
-
]
|
| 260 |
-
},
|
| 261 |
-
{
|
| 262 |
-
"cell_type": "code",
|
| 263 |
-
"execution_count": null,
|
| 264 |
-
"id": "1c318f39-7ca0-4f3c-80fb-73f72202e4e0",
|
| 265 |
-
"metadata": {},
|
| 266 |
-
"outputs": [],
|
| 267 |
-
"source": [
|
| 268 |
-
"%%time\n",
|
| 269 |
-
"tif_file = 'natcrop_fii_100m_cog.tif'\n",
|
| 270 |
-
"vec_file = './cpad-stats-temp.parquet'\n",
|
| 271 |
-
"df = big_zonal_stats(vec_file, tif_file, stats = ['mean'],\n",
|
| 272 |
-
" col_name = \"forest_integrity_loss\", n_jobs=threads, verbose=0).to_parquet(\"cpad-stats-temp.parquet\")\n",
|
| 273 |
-
"\n"
|
| 274 |
-
]
|
| 275 |
-
},
|
| 276 |
-
{
|
| 277 |
-
"cell_type": "code",
|
| 278 |
-
"execution_count": null,
|
| 279 |
-
"id": "aef9070a-c87a-463e-81b8-3cc6c5c9d484",
|
| 280 |
-
"metadata": {},
|
| 281 |
-
"outputs": [],
|
| 282 |
-
"source": [
|
| 283 |
-
"%%time\n",
|
| 284 |
-
"tif_file = 'natcrop_expansion_100m_cog.tif'\n",
|
| 285 |
-
"vec_file = './cpad-stats-temp.parquet'\n",
|
| 286 |
-
"df = big_zonal_stats(vec_file, tif_file, stats = ['mean'], col_name = \"crop_expansion\", n_jobs=threads, verbose=0)\n",
|
| 287 |
-
"gpd.GeoDataFrame(df, geometry=\"geometry\").to_parquet(\"cpad-stats-temp.parquet\")\n"
|
| 288 |
-
]
|
| 289 |
-
},
|
| 290 |
-
{
|
| 291 |
-
"cell_type": "code",
|
| 292 |
-
"execution_count": null,
|
| 293 |
-
"id": "d94f937b-b32c-4de1-b4ac-93ce33f0919f",
|
| 294 |
-
"metadata": {},
|
| 295 |
-
"outputs": [],
|
| 296 |
-
"source": [
|
| 297 |
-
"%%time\n",
|
| 298 |
-
"tif_file = 'natcrop_reduction_100m_cog.tif'\n",
|
| 299 |
-
"vec_file = './cpad-stats-temp.parquet'\n",
|
| 300 |
-
"df = big_zonal_stats(vec_file, tif_file, stats = ['mean'], col_name = \"crop_reduction\", n_jobs=threads, verbose=0).to_parquet(\"cpad-stats-temp.parquet\")\n"
|
| 301 |
-
]
|
| 302 |
-
},
|
| 303 |
-
{
|
| 304 |
-
"cell_type": "code",
|
| 305 |
-
"execution_count": null,
|
| 306 |
-
"id": "6bdaba61-30c1-49d6-a4e6-db68f1daafa3",
|
| 307 |
-
"metadata": {},
|
| 308 |
-
"outputs": [],
|
| 309 |
-
"source": [
|
| 310 |
-
"%%time\n",
|
| 311 |
-
"tif_file = 'hfp_2021_100m_v1-2_cog.tif'\n",
|
| 312 |
-
"vec_file = './cpad-stats-temp.parquet'\n",
|
| 313 |
-
"df = big_zonal_stats(vec_file, tif_file, stats = ['mean'], col_name = \"human_impact\", n_jobs=threads, verbose=0).to_parquet(\"cpad-stats-temp.parquet\")\n"
|
| 314 |
-
]
|
| 315 |
-
},
|
| 316 |
-
{
|
| 317 |
-
"cell_type": "markdown",
|
| 318 |
-
"id": "f8e037d4-7a34-42bc-941f-0c09ee80ef3b",
|
| 319 |
-
"metadata": {},
|
| 320 |
-
"source": [
|
| 321 |
-
"# Need to convert SVI & Justice40 files to tif"
|
| 322 |
-
]
|
| 323 |
-
},
|
| 324 |
-
{
|
| 325 |
-
"cell_type": "code",
|
| 326 |
-
"execution_count": null,
|
| 327 |
-
"id": "c4a19013-65f1-4eef-be2d-0cf1be3d0f7f",
|
| 328 |
-
"metadata": {},
|
| 329 |
-
"outputs": [],
|
| 330 |
-
"source": [
|
| 331 |
-
"import geopandas as gpd\n",
|
| 332 |
-
"import numpy as np\n",
|
| 333 |
-
"import rasterio\n",
|
| 334 |
-
"from rasterio.features import rasterize\n",
|
| 335 |
-
"from rasterio.transform import from_bounds\n",
|
| 336 |
-
"\n",
|
| 337 |
-
"def get_geotiff(gdf, output_file,col):\n",
|
| 338 |
-
" gdf = gdf.set_geometry(\"geometry\")\n",
|
| 339 |
-
" gdf = gdf.set_crs(\"EPSG:4326\")\n",
|
| 340 |
-
" print(gdf.crs)\n",
|
| 341 |
-
"\n",
|
| 342 |
-
" # Set raster properties\n",
|
| 343 |
-
" minx, miny, maxx, maxy = gdf.total_bounds # Get the bounds of the geometry\n",
|
| 344 |
-
" pixel_size = 0.01 # Define the pixel size in units of the CRS\n",
|
| 345 |
-
" width = int((maxx - minx) / pixel_size)\n",
|
| 346 |
-
" height = int((maxy - miny) / pixel_size)\n",
|
| 347 |
-
" transform = from_bounds(minx, miny, maxx, maxy, width, height)\n",
|
| 348 |
-
" \n",
|
| 349 |
-
" # Define rasterization with continuous values\n",
|
| 350 |
-
" shapes = ((geom, value) for geom, value in zip(gdf.geometry, gdf[col]))\n",
|
| 351 |
-
" raster = rasterize(\n",
|
| 352 |
-
" shapes,\n",
|
| 353 |
-
" out_shape=(height, width),\n",
|
| 354 |
-
" transform=transform,\n",
|
| 355 |
-
" fill=0.0, # Background value for areas outside the geometry\n",
|
| 356 |
-
" dtype=\"float32\" # Set data type to handle continuous values\n",
|
| 357 |
-
" )\n",
|
| 358 |
-
" print(\"Unique values in raster:\", np.unique(raster))\n",
|
| 359 |
-
"\n",
|
| 360 |
-
" # Define GeoTIFF metadata\n",
|
| 361 |
-
" out_meta = {\n",
|
| 362 |
-
" \"driver\": \"GTiff\",\n",
|
| 363 |
-
" \"height\": height,\n",
|
| 364 |
-
" \"width\": width,\n",
|
| 365 |
-
" \"count\": 1,\n",
|
| 366 |
-
" \"dtype\": raster.dtype,\n",
|
| 367 |
-
" \"crs\": gdf.crs,\n",
|
| 368 |
-
" \"transform\": transform,\n",
|
| 369 |
-
" \"compress\": \"deflate\" # Use compression to reduce file size\n",
|
| 370 |
-
" }\n",
|
| 371 |
-
" \n",
|
| 372 |
-
" # Write to a GeoTIFF file with COG options\n",
|
| 373 |
-
" with rasterio.open(output_file, \"w\", **out_meta) as dest:\n",
|
| 374 |
-
" dest.write(raster, 1)\n",
|
| 375 |
-
" dest.build_overviews([2, 4, 8, 16], rasterio.enums.Resampling.average)\n",
|
| 376 |
-
" dest.update_tags(1, TIFFTAG_RESOLUTION_UNIT=\"Meter\")\n"
|
| 377 |
-
]
|
| 378 |
-
},
|
| 379 |
-
{
|
| 380 |
-
"cell_type": "markdown",
|
| 381 |
-
"id": "f4925a74-5ed2-49a4-845b-6a0f0398a43e",
|
| 382 |
-
"metadata": {},
|
| 383 |
-
"source": [
|
| 384 |
-
"# SVI"
|
| 385 |
-
]
|
| 386 |
-
},
|
| 387 |
-
{
|
| 388 |
-
"cell_type": "code",
|
| 389 |
-
"execution_count": null,
|
| 390 |
-
"id": "4e678f01-73af-4f99-a565-e9b7f04d9547",
|
| 391 |
-
"metadata": {},
|
| 392 |
-
"outputs": [],
|
| 393 |
-
"source": [
|
| 394 |
-
"# clean up SVI data\n",
|
| 395 |
-
"svi_df = (con\n",
|
| 396 |
-
" .read_parquet(\"svi2020_us_tract.parquet\")\n",
|
| 397 |
-
" .select(\"RPL_THEMES\",\"RPL_THEME1\",\"RPL_THEME2\",\"RPL_THEME3\",\"RPL_THEME4\",\"Shape\")\n",
|
| 398 |
-
" .rename(SVI = \"RPL_THEMES\", socioeconomic = \"RPL_THEME1\", \n",
|
| 399 |
-
" household_char = \"RPL_THEME2\", racial_ethnic_minority = \"RPL_THEME3\",\n",
|
| 400 |
-
" housing_transit = \"RPL_THEME4\", geometry = \"Shape\")\n",
|
| 401 |
-
".cast({\"geometry\":\"geometry\"})\n",
|
| 402 |
-
")\n",
|
| 403 |
-
"svi_df.execute().to_parquet(\"svi2020_us_tract_clean.parquet\")\n"
|
| 404 |
-
]
|
| 405 |
-
},
|
| 406 |
-
{
|
| 407 |
-
"cell_type": "code",
|
| 408 |
-
"execution_count": null,
|
| 409 |
-
"id": "c5046d6b-9798-46d3-a1bc-548e29414007",
|
| 410 |
-
"metadata": {},
|
| 411 |
-
"outputs": [],
|
| 412 |
-
"source": [
|
| 413 |
-
"gdf = gpd.read_parquet(\"svi2020_us_tract_clean.parquet\")\n",
|
| 414 |
-
"svi = gdf[['SVI','geometry']]\n",
|
| 415 |
-
"socio = gdf[['socioeconomic','geometry']]\n",
|
| 416 |
-
"house = gdf[['household_char','geometry']]\n",
|
| 417 |
-
"minority = gdf[['racial_ethnic_minority','geometry']]\n",
|
| 418 |
-
"transit = gdf[['housing_transit','geometry']]\n",
|
| 419 |
-
"\n",
|
| 420 |
-
"#convert SVI parquet to tif\n",
|
| 421 |
-
"get_geotiff(svi,\"svi.tif\",\"SVI\")\n",
|
| 422 |
-
"get_geotiff(socio,\"svi_socioeconomic.tif\",\"socioeconomic\")\n",
|
| 423 |
-
"get_geotiff(house,\"svi_household.tif\",\"household_char\")\n",
|
| 424 |
-
"get_geotiff(minority,\"svi_minority.tif\",\"racial_ethnic_minority\")\n",
|
| 425 |
-
"get_geotiff(transit,\"svi_transit.tif\",\"housing_transit\")"
|
| 426 |
-
]
|
| 427 |
-
},
|
| 428 |
-
{
|
| 429 |
-
"cell_type": "code",
|
| 430 |
-
"execution_count": null,
|
| 431 |
-
"id": "6a36b77f-d0be-45bd-9318-da4b57eaf353",
|
| 432 |
-
"metadata": {},
|
| 433 |
-
"outputs": [],
|
| 434 |
-
"source": [
|
| 435 |
-
"%%time\n",
|
| 436 |
-
"tif_file = 'svi.tif'\n",
|
| 437 |
-
"vec_file = './cpad-stats-temp.parquet'\n",
|
| 438 |
-
"df = big_zonal_stats(vec_file, tif_file, stats = ['mean'], col_name = \"SVI\", n_jobs=threads, verbose=0).to_parquet(\"cpad-stats-temp.parquet\")\n",
|
| 439 |
-
"\n"
|
| 440 |
-
]
|
| 441 |
-
},
|
| 442 |
-
{
|
| 443 |
-
"cell_type": "code",
|
| 444 |
-
"execution_count": null,
|
| 445 |
-
"id": "05ef74e2-3f23-4f69-8cd3-8862cb73a259",
|
| 446 |
-
"metadata": {},
|
| 447 |
-
"outputs": [],
|
| 448 |
-
"source": [
|
| 449 |
-
"%%time\n",
|
| 450 |
-
"vec_file = './cpad-stats-temp.parquet'\n",
|
| 451 |
-
"tif_file = 'svi_socioeconomic.tif'\n",
|
| 452 |
-
"df = big_zonal_stats(vec_file, tif_file, stats = ['mean'], col_name = \"socioeconomic_status\", n_jobs=threads, verbose=0).to_parquet(\"cpad-stats-temp.parquet\")\n",
|
| 453 |
-
"\n"
|
| 454 |
-
]
|
| 455 |
-
},
|
| 456 |
-
{
|
| 457 |
-
"cell_type": "code",
|
| 458 |
-
"execution_count": null,
|
| 459 |
-
"id": "23417a03-38c2-4b31-8340-f08ec8a11631",
|
| 460 |
-
"metadata": {},
|
| 461 |
-
"outputs": [],
|
| 462 |
-
"source": [
|
| 463 |
-
"%%time\n",
|
| 464 |
-
"vec_file = './cpad-stats-temp.parquet'\n",
|
| 465 |
-
"tif_file = 'svi_household.tif'\n",
|
| 466 |
-
"df = big_zonal_stats(vec_file, tif_file, stats = ['mean'], col_name = \"household_char\", n_jobs=threads, verbose=0).to_parquet(\"cpad-stats-temp.parquet\")\n",
|
| 467 |
-
"\n"
|
| 468 |
-
]
|
| 469 |
-
},
|
| 470 |
-
{
|
| 471 |
-
"cell_type": "code",
|
| 472 |
-
"execution_count": null,
|
| 473 |
-
"id": "de86d7f0-6cdc-4d05-bdee-d9803cbd83bd",
|
| 474 |
-
"metadata": {},
|
| 475 |
-
"outputs": [],
|
| 476 |
-
"source": [
|
| 477 |
-
"%%time\n",
|
| 478 |
-
"vec_file = './cpad-stats-temp.parquet'\n",
|
| 479 |
-
"tif_file = 'svi_minority.tif'\n",
|
| 480 |
-
"df = big_zonal_stats(vec_file, tif_file, stats = ['mean'], col_name = \"racial_ethnic_minority\", n_jobs=threads, verbose=0).to_parquet(\"cpad-stats-temp.parquet\")\n"
|
| 481 |
-
]
|
| 482 |
-
},
|
| 483 |
-
{
|
| 484 |
-
"cell_type": "code",
|
| 485 |
-
"execution_count": null,
|
| 486 |
-
"id": "0c49dd50-7dd3-4240-9af8-3e32ec656bc0",
|
| 487 |
-
"metadata": {},
|
| 488 |
-
"outputs": [],
|
| 489 |
-
"source": [
|
| 490 |
-
"%%time\n",
|
| 491 |
-
"vec_file = './cpad-stats-temp.parquet'\n",
|
| 492 |
-
"tif_file = 'svi_transit.tif'\n",
|
| 493 |
-
"df = big_zonal_stats(vec_file, tif_file, stats = ['mean'], col_name = \"housing_transit\", n_jobs=threads, verbose=0).to_parquet(\"cpad-stats-temp.parquet\")\n"
|
| 494 |
-
]
|
| 495 |
-
},
|
| 496 |
-
{
|
| 497 |
-
"cell_type": "markdown",
|
| 498 |
-
"id": "ff4b6604-9828-4882-90bd-554c21f5c6e6",
|
| 499 |
-
"metadata": {},
|
| 500 |
-
"source": [
|
| 501 |
-
"# Justice40 "
|
| 502 |
-
]
|
| 503 |
-
},
|
| 504 |
-
{
|
| 505 |
-
"cell_type": "code",
|
| 506 |
-
"execution_count": null,
|
| 507 |
-
"id": "3678a91f-72f7-4339-a409-a97776cba043",
|
| 508 |
-
"metadata": {},
|
| 509 |
-
"outputs": [],
|
| 510 |
-
"source": [
|
| 511 |
-
"#clean up\n",
|
| 512 |
-
"justice40 = (con\n",
|
| 513 |
-
" .read_parquet(\"disadvantaged-communities.parquet\")\n",
|
| 514 |
-
" .rename(geometry = \"SHAPE\",justice40=\"Disadvan\")\n",
|
| 515 |
-
" .filter(_.StateName == \"California\")\n",
|
| 516 |
-
" .mutate(geometry = _.geometry.convert(\"ESRI:102039\",\"EPSG:4326\"))\n",
|
| 517 |
-
" .select(\"justice40\",\"geometry\")\n",
|
| 518 |
-
" )\n",
|
| 519 |
-
"justice40.execute().to_parquet(\"ca_justice40.parquet\")"
|
| 520 |
-
]
|
| 521 |
-
},
|
| 522 |
-
{
|
| 523 |
-
"cell_type": "code",
|
| 524 |
-
"execution_count": null,
|
| 525 |
-
"id": "8faa425f-6f9c-4189-a53a-24dd0250c539",
|
| 526 |
-
"metadata": {},
|
| 527 |
-
"outputs": [],
|
| 528 |
-
"source": [
|
| 529 |
-
"# #justice40 is either 0 or 1, so we want to get the percentage of polygon where justice40 = 1. \n",
|
| 530 |
-
"\n",
|
| 531 |
-
"def big_zonal_stats_binary(vec_file, justice40_file, col_name,projected_crs=\"EPSG:3310\"):\n",
|
| 532 |
-
" # Read both vector files as GeoDataFrames\n",
|
| 533 |
-
" gdf = gpd.read_parquet(vec_file)\n",
|
| 534 |
-
" justice40_gdf = gpd.read_parquet(justice40_file)\n",
|
| 535 |
-
" \n",
|
| 536 |
-
" # Set CRS if not already set (assuming both should be in EPSG:4326, modify if needed)\n",
|
| 537 |
-
" if gdf.crs is None:\n",
|
| 538 |
-
" gdf = gdf.set_crs(\"EPSG:4326\")\n",
|
| 539 |
-
" if justice40_gdf.crs is None:\n",
|
| 540 |
-
" justice40_gdf = justice40_gdf.set_crs(\"EPSG:4326\")\n",
|
| 541 |
-
" # Ensure both GeoDataFrames are in the same CRS and reproject to a projected CRS for area calculations\n",
|
| 542 |
-
" gdf = gdf.to_crs(projected_crs)\n",
|
| 543 |
-
" justice40_gdf = justice40_gdf.to_crs(projected_crs)\n",
|
| 544 |
-
" \n",
|
| 545 |
-
" # Ensure both GeoDataFrames are in the same CRS\n",
|
| 546 |
-
" gdf = gdf.to_crs(justice40_gdf.crs)\n",
|
| 547 |
-
" \n",
|
| 548 |
-
" # Filter justice40 polygons where justice40 == 1\n",
|
| 549 |
-
" justice40_gdf = justice40_gdf[justice40_gdf['justice40'] == 1].copy()\n",
|
| 550 |
-
" \n",
|
| 551 |
-
" # Prepare a list to hold percentage of justice40 == 1 for each polygon\n",
|
| 552 |
-
" percentages = []\n",
|
| 553 |
-
" \n",
|
| 554 |
-
" # Iterate over each polygon in the main GeoDataFrame\n",
|
| 555 |
-
" for geom in gdf.geometry:\n",
|
| 556 |
-
" # Find intersecting justice40 polygons\n",
|
| 557 |
-
" justice40_intersections = justice40_gdf[justice40_gdf.intersects(geom)].copy()\n",
|
| 558 |
-
" \n",
|
| 559 |
-
" # Calculate the intersection area\n",
|
| 560 |
-
" if not justice40_intersections.empty:\n",
|
| 561 |
-
" justice40_intersections['intersection'] = justice40_intersections.intersection(geom)\n",
|
| 562 |
-
" total_intersection_area = justice40_intersections['intersection'].area.sum()\n",
|
| 563 |
-
" \n",
|
| 564 |
-
" # Calculate percentage based on original polygon's area\n",
|
| 565 |
-
" percentage_1 = (total_intersection_area / geom.area) \n",
|
| 566 |
-
" else:\n",
|
| 567 |
-
" percentage_1 = 0.0 # No intersection with justice40 == 1 polygons\n",
|
| 568 |
-
" \n",
|
| 569 |
-
" # Append result\n",
|
| 570 |
-
" percentages.append(percentage_1)\n",
|
| 571 |
-
" \n",
|
| 572 |
-
" # Add results to the original GeoDataFrame\n",
|
| 573 |
-
" gdf[col_name] = percentages\n",
|
| 574 |
-
" return gdf\n",
|
| 575 |
-
"\n",
|
| 576 |
-
"\n"
|
| 577 |
-
]
|
| 578 |
-
},
|
| 579 |
-
{
|
| 580 |
-
"cell_type": "code",
|
| 581 |
-
"execution_count": null,
|
| 582 |
-
"id": "fe80fc28-73ce-4a26-9925-851c2798e467",
|
| 583 |
-
"metadata": {},
|
| 584 |
-
"outputs": [],
|
| 585 |
-
"source": [
|
| 586 |
-
"%%time\n",
|
| 587 |
-
"vec_file = './cpad-stats-temp.parquet'\n",
|
| 588 |
-
"\n",
|
| 589 |
-
"df = big_zonal_stats_binary(vec_file, \"ca_justice40.parquet\", col_name=\"percent_disadvantaged\")\n",
|
| 590 |
-
"df.to_parquet(\"cpad-stats-temp.parquet\")\n"
|
| 591 |
-
]
|
| 592 |
-
},
|
| 593 |
-
{
|
| 594 |
-
"cell_type": "markdown",
|
| 595 |
-
"id": "5438a4f4-377e-41fe-800b-8ffc1f33caa0",
|
| 596 |
-
"metadata": {},
|
| 597 |
-
"source": [
|
| 598 |
-
"# Fire"
|
| 599 |
-
]
|
| 600 |
-
},
|
| 601 |
-
{
|
| 602 |
-
"cell_type": "code",
|
| 603 |
-
"execution_count": null,
|
| 604 |
-
"id": "4bd83b4d-01df-49d8-99e1-6740d365c833",
|
| 605 |
-
"metadata": {},
|
| 606 |
-
"outputs": [],
|
| 607 |
-
"source": [
|
| 608 |
-
"import geopandas as gpd\n",
|
| 609 |
-
"\n",
|
| 610 |
-
"#get percentage of polygon with fire occurrence \n",
|
| 611 |
-
"def fire_stats(file_name, fire_df, col_name):\n",
|
| 612 |
-
" gdf = gpd.read_parquet(file_name)\n",
|
| 613 |
-
" \n",
|
| 614 |
-
" percentages = []\n",
|
| 615 |
-
" # Find all fires that intersect with the current protected area \n",
|
| 616 |
-
" for geom in gdf.geometry:\n",
|
| 617 |
-
" fire_intersections = fire_df[fire_df.intersects(geom)].copy()\n",
|
| 618 |
-
" if not fire_intersections.empty:\n",
|
| 619 |
-
" # If there is only one intersecting fire, compute the intersection area\n",
|
| 620 |
-
" if len(fire_intersections) == 1:\n",
|
| 621 |
-
" intersection_area = fire_intersections.geometry.iloc[0].intersection(geom).area\n",
|
| 622 |
-
" else:\n",
|
| 623 |
-
" # If there are multiple intersecting fires, use a union to avoid double-counting\n",
|
| 624 |
-
" unioned_fires = fire_intersections.unary_union\n",
|
| 625 |
-
" intersection_area = unioned_fires.intersection(geom).area\n",
|
| 626 |
-
" \n",
|
| 627 |
-
" percentage_1 = round((intersection_area / geom.area),3)\n",
|
| 628 |
-
" else:\n",
|
| 629 |
-
" percentage_1 = 0.0 \n",
|
| 630 |
-
"\n",
|
| 631 |
-
" percentages.append(percentage_1)\n",
|
| 632 |
-
" \n",
|
| 633 |
-
" gdf[col_name] = percentages\n",
|
| 634 |
-
" return gdf\n"
|
| 635 |
-
]
|
| 636 |
-
},
|
| 637 |
-
{
|
| 638 |
-
"cell_type": "code",
|
| 639 |
-
"execution_count": null,
|
| 640 |
-
"id": "4ce35cea-8897-42c0-b1f6-01b414a5b556",
|
| 641 |
-
"metadata": {},
|
| 642 |
-
"outputs": [],
|
| 643 |
-
"source": [
|
| 644 |
-
"#historical fire perimeters \n",
|
| 645 |
-
"fire_20 = (con\n",
|
| 646 |
-
" .read_parquet(\"firep22_1.parquet\")\n",
|
| 647 |
-
" .rename(year = \"YEAR_\")\n",
|
| 648 |
-
" .filter(_.STATE == \"CA\", _.year != '')\n",
|
| 649 |
-
" .cast({\"year\":\"int\"})\n",
|
| 650 |
-
" .filter(_.year>=2003)\n",
|
| 651 |
-
" .select(\"year\",\"geometry\")\n",
|
| 652 |
-
" .mutate(\n",
|
| 653 |
-
" geometry=ibis.ifelse(\n",
|
| 654 |
-
" _.geometry.is_valid(),\n",
|
| 655 |
-
" _.geometry, # Keep the geometry if it's valid\n",
|
| 656 |
-
" _.geometry.buffer(0) # Apply buffer(0) to fix invalid geometries\n",
|
| 657 |
-
" )\n",
|
| 658 |
-
" )\n",
|
| 659 |
-
" )\n",
|
| 660 |
-
"fire_20.execute().to_parquet(\"ca-fire-20yrs.parquet\")\n",
|
| 661 |
-
"fire_10 = fire_20.filter(_.year>=2013)\n",
|
| 662 |
-
"fire_5 = fire_20.filter(_.year>=2018)\n",
|
| 663 |
-
"fire_2 = fire_20.filter(_.year>=2022)\n",
|
| 664 |
-
"\n",
|
| 665 |
-
"\n",
|
| 666 |
-
"fire_20_df = fire_20.execute().set_crs(\"EPSG:3310\")\n",
|
| 667 |
-
"fire_10_df = fire_10.execute().set_crs(\"EPSG:3310\")\n",
|
| 668 |
-
"fire_5_df = fire_5.execute().set_crs(\"EPSG:3310\")\n",
|
| 669 |
-
"fire_2_df = fire_2.execute().set_crs(\"EPSG:3310\")\n"
|
| 670 |
-
]
|
| 671 |
-
},
|
| 672 |
-
{
|
| 673 |
-
"cell_type": "code",
|
| 674 |
-
"execution_count": null,
|
| 675 |
-
"id": "0a041210-6ffe-49b0-b4a7-3a9220acedb9",
|
| 676 |
-
"metadata": {},
|
| 677 |
-
"outputs": [],
|
| 678 |
-
"source": [
|
| 679 |
-
"#prescribed burns\n",
|
| 680 |
-
"rxburn_20 = (con\n",
|
| 681 |
-
" .read_parquet(\"rxburn22_1.parquet\")\n",
|
| 682 |
-
" .rename(year = \"YEAR_\")\n",
|
| 683 |
-
" .filter(_.STATE == \"CA\", _.year != ' ', _.year != '')\n",
|
| 684 |
-
" .cast({\"year\":\"int\"})\n",
|
| 685 |
-
" .filter(_.year>=2003)\n",
|
| 686 |
-
" .select(\"year\",\"geometry\")\n",
|
| 687 |
-
" .mutate(\n",
|
| 688 |
-
" geometry=ibis.ifelse(\n",
|
| 689 |
-
" _.geometry.is_valid(),\n",
|
| 690 |
-
" _.geometry, # Keep the geometry if it's valid\n",
|
| 691 |
-
" _.geometry.buffer(0) # Apply buffer(0) to fix invalid geometries\n",
|
| 692 |
-
" )\n",
|
| 693 |
-
" )\n",
|
| 694 |
-
" )\n",
|
| 695 |
-
"\n",
|
| 696 |
-
"rxburn_20.execute().to_parquet(\"ca-rxburn-20yrs.parquet\")\n",
|
| 697 |
-
"rxburn_10 = (rxburn_20.filter(_.year>=2013))\n",
|
| 698 |
-
"rxburn_5 = (rxburn_20.filter(_.year>=2018))\n",
|
| 699 |
-
"rxburn_2 = (rxburn_20.filter(_.year>=2022))\n",
|
| 700 |
-
"\n",
|
| 701 |
-
"rxburn_20_df = rxburn_20.execute().set_crs(\"EPSG:3310\")\n",
|
| 702 |
-
"rxburn_10_df = rxburn_10.execute().set_crs(\"EPSG:3310\")\n",
|
| 703 |
-
"rxburn_5_df = rxburn_5.execute().set_crs(\"EPSG:3310\")\n",
|
| 704 |
-
"rxburn_2_df = rxburn_2.execute().set_crs(\"EPSG:3310\")"
|
| 705 |
-
]
|
| 706 |
-
},
|
| 707 |
-
{
|
| 708 |
-
"cell_type": "code",
|
| 709 |
-
"execution_count": null,
|
| 710 |
-
"id": "fc955b02-efc1-4ae3-b8e4-ea424d491a68",
|
| 711 |
-
"metadata": {},
|
| 712 |
-
"outputs": [],
|
| 713 |
-
"source": [
|
| 714 |
-
"# need to validate geometries, using epsg:3310 to match fire polygons\n",
|
| 715 |
-
"ca = (con\n",
|
| 716 |
-
" .read_parquet('cpad-stats-temp.parquet')\n",
|
| 717 |
-
" .mutate(geom = _.geom.convert(\"EPSG:4326\",\"EPSG:3310\"))\n",
|
| 718 |
-
" .mutate(\n",
|
| 719 |
-
" geometry=ibis.ifelse(\n",
|
| 720 |
-
" _.geom.is_valid(),\n",
|
| 721 |
-
" _.geom, # Keep the geometry if it's valid\n",
|
| 722 |
-
" _.geom.buffer(0) # Apply buffer(0) to fix invalid geometries\n",
|
| 723 |
-
" )\n",
|
| 724 |
-
" )\n",
|
| 725 |
-
" .drop('geom')\n",
|
| 726 |
-
" )\n",
|
| 727 |
-
"gdf = ca.execute()\n",
|
| 728 |
-
"gdf = gdf.set_crs('EPSG:3310')\n",
|
| 729 |
-
"gdf.to_parquet('cpad-stats-temp-EPSG3310.parquet')\n"
|
| 730 |
-
]
|
| 731 |
-
},
|
| 732 |
-
{
|
| 733 |
-
"cell_type": "code",
|
| 734 |
-
"execution_count": null,
|
| 735 |
-
"id": "68e25266-efc8-4378-afc5-95c7a769ca81",
|
| 736 |
-
"metadata": {},
|
| 737 |
-
"outputs": [],
|
| 738 |
-
"source": [
|
| 739 |
-
"%%time\n",
|
| 740 |
-
"file_name = 'cpad-stats-temp-EPSG3310.parquet'\n",
|
| 741 |
-
"\n",
|
| 742 |
-
"names = [\"percent_fire_20yr\", \"percent_fire_10yr\", \"percent_fire_5yr\",\n",
|
| 743 |
-
" \"percent_fire_2yr\",\"percent_rxburn_20yr\", \"percent_rxburn_10yr\", \n",
|
| 744 |
-
" \"percent_rxburn_5yr\",\"percent_rxburn_2yr\"]\n",
|
| 745 |
-
"dfs = [fire_20_df,fire_10_df,fire_5_df,fire_2_df,rxburn_20_df,rxburn_10_df,rxburn_5_df,rxburn_2_df]\n",
|
| 746 |
-
"\n",
|
| 747 |
-
"for df,name in zip(dfs,names):\n",
|
| 748 |
-
" df_stat = fire_stats(file_name,df, col_name=name)\n",
|
| 749 |
-
" df_stat.to_parquet(file_name)"
|
| 750 |
-
]
|
| 751 |
-
},
|
| 752 |
-
{
|
| 753 |
-
"cell_type": "code",
|
| 754 |
-
"execution_count": null,
|
| 755 |
-
"id": "cd4acb35-d1a3-4632-ae30-c6e3e923e94c",
|
| 756 |
-
"metadata": {},
|
| 757 |
-
"outputs": [],
|
| 758 |
-
"source": [
|
| 759 |
-
"#save data back to cpad-stats-temp\n",
|
| 760 |
-
"# (not really necessary but I want to reuse the same code)\n",
|
| 761 |
-
"ca = (con\n",
|
| 762 |
-
" .read_parquet(file_name)\n",
|
| 763 |
-
" .mutate(geometry = _.geometry.convert(\"EPSG:3310\",\"EPSG:4326\"))\n",
|
| 764 |
-
" )\n",
|
| 765 |
-
"gdf = ca.execute()\n",
|
| 766 |
-
"gdf= gdf.set_crs('EPSG:4326')\n",
|
| 767 |
-
"gdf.to_parquet(\"cpad-stats-temp.parquet\")\n",
|
| 768 |
-
"\n"
|
| 769 |
-
]
|
| 770 |
-
},
|
| 771 |
-
{
|
| 772 |
-
"cell_type": "markdown",
|
| 773 |
-
"id": "e3083b85-1322-4188-ac08-e73c2570978c",
|
| 774 |
-
"metadata": {},
|
| 775 |
-
"source": [
|
| 776 |
-
"# Cleaning up + Rounding floats"
|
| 777 |
-
]
|
| 778 |
-
},
|
| 779 |
-
{
|
| 780 |
-
"cell_type": "code",
|
| 781 |
-
"execution_count": null,
|
| 782 |
-
"id": "2e4de199-82d4-4e2b-8572-6fe19b57d1ee",
|
| 783 |
-
"metadata": {},
|
| 784 |
-
"outputs": [],
|
| 785 |
-
"source": [
|
| 786 |
-
"## clean up\n",
|
| 787 |
-
"con = ibis.duckdb.connect(extensions=[\"spatial\"])\n",
|
| 788 |
-
"ca_geom = con.read_parquet(\"ca2024-30m.parquet\").cast({\"geom\":\"geometry\"}).select(\"id\",\"geom\")\n",
|
| 789 |
-
"\n",
|
| 790 |
-
"ca = (con\n",
|
| 791 |
-
" .read_parquet(\"cpad-stats-temp.parquet\")\n",
|
| 792 |
-
" .cast({\n",
|
| 793 |
-
" \"crop_expansion\": \"int64\",\n",
|
| 794 |
-
" \"crop_reduction\": \"int64\",\n",
|
| 795 |
-
" \"manageable_carbon\": \"int64\",\n",
|
| 796 |
-
" \"irrecoverable_carbon\": \"int64\"\n",
|
| 797 |
-
" })\n",
|
| 798 |
-
" .mutate(\n",
|
| 799 |
-
" richness=_.richness.round(3),\n",
|
| 800 |
-
" rsr=_.rsr.round(3),\n",
|
| 801 |
-
" all_species_rwr=_.all_species_rwr.round(3),\n",
|
| 802 |
-
" all_species_richness=_.all_species_richness.round(3),\n",
|
| 803 |
-
" percent_disadvantaged=(_.percent_disadvantaged).round(3),\n",
|
| 804 |
-
" svi=_.svi.round(3),\n",
|
| 805 |
-
" svi_socioeconomic_status=_.socioeconomic_status.round(3),\n",
|
| 806 |
-
" svi_household_char=_.household_char.round(3),\n",
|
| 807 |
-
" svi_racial_ethnic_minority=_.racial_ethnic_minority.round(3),\n",
|
| 808 |
-
" svi_housing_transit=_.housing_transit.round(3),\n",
|
| 809 |
-
" human_impact=_.human_impact.round(3),\n",
|
| 810 |
-
" deforest_carbon=_.deforest_carbon.round(3),\n",
|
| 811 |
-
" biodiversity_intactness_loss=_.biodiversity_intactness_loss.round(3),\n",
|
| 812 |
-
" forest_integrity_loss=_.forest_integrity_loss.round(3),\n",
|
| 813 |
-
" percent_fire_20yr = _.percent_fire_20yr.round(3),\n",
|
| 814 |
-
" percent_fire_10yr = _.percent_fire_10yr.round(3),\n",
|
| 815 |
-
" percent_fire_5yr = _.percent_fire_5yr.round(3),\n",
|
| 816 |
-
" percent_fire_2yr = _.percent_fire_2yr.round(3),\n",
|
| 817 |
-
" percent_rxburn_20yr = _.percent_rxburn_20yr.round(3),\n",
|
| 818 |
-
" percent_rxburn_10yr = _.percent_rxburn_10yr.round(3),\n",
|
| 819 |
-
" percent_rxburn_5yr = _.percent_rxburn_5yr.round(3),\n",
|
| 820 |
-
" percent_rxburn_2yr = _.percent_rxburn_2yr.round(3),\n",
|
| 821 |
-
" )\n",
|
| 822 |
-
" # only grabbing columns we are making charts with \n",
|
| 823 |
-
" .select('established', 'reGAP', 'name', 'access_type', 'manager', 'manager_type', 'Easement', 'Acres', 'id', 'type','richness', \n",
|
| 824 |
-
" 'rsr', 'irrecoverable_carbon', 'manageable_carbon', 'percent_fire_20yr', 'percent_fire_10yr', 'percent_fire_5yr','percent_fire_2yr',\n",
|
| 825 |
-
" 'percent_rxburn_20yr', 'percent_rxburn_10yr', 'percent_rxburn_5yr','percent_rxburn_2yr', 'percent_disadvantaged',\n",
|
| 826 |
-
" 'svi', 'svi_socioeconomic_status', 'svi_household_char', 'svi_racial_ethnic_minority',\n",
|
| 827 |
-
" 'svi_housing_transit', 'deforest_carbon','human_impact'\n",
|
| 828 |
-
" )\n",
|
| 829 |
-
" .join(ca_geom, \"id\", how=\"inner\")\n",
|
| 830 |
-
" )\n",
|
| 831 |
-
"\n",
|
| 832 |
-
"ca.head(5).execute()\n"
|
| 833 |
-
]
|
| 834 |
-
},
|
| 835 |
-
{
|
| 836 |
-
"cell_type": "markdown",
|
| 837 |
-
"id": "3780de2c-3a68-442c-bb3b-64c792418979",
|
| 838 |
-
"metadata": {},
|
| 839 |
-
"source": [
|
| 840 |
-
"# Save as PMTiles + Upload data"
|
| 841 |
-
]
|
| 842 |
-
},
|
| 843 |
-
{
|
| 844 |
-
"cell_type": "code",
|
| 845 |
-
"execution_count": 1,
|
| 846 |
-
"id": "05c791c9-888a-483a-9dbb-a2ba7eb1bce2",
|
| 847 |
-
"metadata": {},
|
| 848 |
-
"outputs": [
|
| 849 |
-
{
|
| 850 |
-
"name": "stderr",
|
| 851 |
-
"output_type": "stream",
|
| 852 |
-
"text": [
|
| 853 |
-
"Note: Environment variable`HF_TOKEN` is set and is the current active token independently from the token you've just configured.\n"
|
| 854 |
-
]
|
| 855 |
-
}
|
| 856 |
-
],
|
| 857 |
-
"source": [
|
| 858 |
-
"import subprocess\n",
|
| 859 |
-
"import os\n",
|
| 860 |
-
"from huggingface_hub import HfApi, login\n",
|
| 861 |
-
"import streamlit as st\n",
|
| 862 |
-
"\n",
|
| 863 |
-
"login(st.secrets[\"HF_TOKEN\"])\n",
|
| 864 |
-
"# api = HfApi(add_to_git_credential=False)\n",
|
| 865 |
-
"api = HfApi()\n",
|
| 866 |
-
"\n",
|
| 867 |
-
"def hf_upload(file, repo_id,repo_type):\n",
|
| 868 |
-
" info = api.upload_file(\n",
|
| 869 |
-
" path_or_fileobj=file,\n",
|
| 870 |
-
" path_in_repo=file,\n",
|
| 871 |
-
" repo_id=repo_id,\n",
|
| 872 |
-
" repo_type=repo_type,\n",
|
| 873 |
-
" )\n",
|
| 874 |
-
"def generate_pmtiles(input_file, output_file, max_zoom=12):\n",
|
| 875 |
-
" # Ensure Tippecanoe is installed\n",
|
| 876 |
-
" if subprocess.call([\"which\", \"tippecanoe\"], stdout=subprocess.DEVNULL) != 0:\n",
|
| 877 |
-
" raise RuntimeError(\"Tippecanoe is not installed or not in PATH\")\n",
|
| 878 |
-
"\n",
|
| 879 |
-
" # Construct the Tippecanoe command\n",
|
| 880 |
-
" command = [\n",
|
| 881 |
-
" \"tippecanoe\",\n",
|
| 882 |
-
" \"-o\", output_file,\n",
|
| 883 |
-
" \"-zg\",\n",
|
| 884 |
-
" \"--extend-zooms-if-still-dropping\",\n",
|
| 885 |
-
" \"--force\",\n",
|
| 886 |
-
" \"--projection\", \"EPSG:4326\", \n",
|
| 887 |
-
" \"-L\",\"layer:\"+input_file,\n",
|
| 888 |
-
" ]\n",
|
| 889 |
-
" # Run Tippecanoe\n",
|
| 890 |
-
" try:\n",
|
| 891 |
-
" subprocess.run(command, check=True)\n",
|
| 892 |
-
" print(f\"Successfully generated PMTiles file: {output_file}\")\n",
|
| 893 |
-
" except subprocess.CalledProcessError as e:\n",
|
| 894 |
-
" print(f\"Error running Tippecanoe: {e}\")\n",
|
| 895 |
-
"\n"
|
| 896 |
-
]
|
| 897 |
-
},
|
| 898 |
-
{
|
| 899 |
-
"cell_type": "code",
|
| 900 |
-
"execution_count": null,
|
| 901 |
-
"id": "1f2d179d-6d47-4e84-83c6-7cb3d969fc00",
|
| 902 |
-
"metadata": {},
|
| 903 |
-
"outputs": [],
|
| 904 |
-
"source": [
|
| 905 |
-
"gdf = ca.execute().set_crs(\"EPSG:4326\")\n",
|
| 906 |
-
"gdf.to_file(\"cpad-stats.geojson\")\n",
|
| 907 |
-
"\n",
|
| 908 |
-
"generate_pmtiles(\"cpad-stats.geojson\", \"cpad-stats.pmtiles\")\n",
|
| 909 |
-
"hf_upload(\"cpad-stats.pmtiles\", \"boettiger-lab/ca-30x30\",\"dataset\")\n",
|
| 910 |
-
"\n",
|
| 911 |
-
"gdf.to_parquet(\"cpad-stats.parquet\")\n",
|
| 912 |
-
"hf_upload(\"cpad-stats.parquet\", \"boettiger-lab/ca-30x30\",\"dataset\")\n",
|
| 913 |
-
"hf_upload(\"cpad-stats.parquet\", \"boettiger-lab/ca-30x30\",\"space\")\n",
|
| 914 |
-
"\n"
|
| 915 |
-
]
|
| 916 |
-
},
|
| 917 |
-
{
|
| 918 |
-
"cell_type": "markdown",
|
| 919 |
-
"id": "09467342-c160-413b-9cdc-31a4bec968cf",
|
| 920 |
-
"metadata": {},
|
| 921 |
-
"source": [
|
| 922 |
-
"# Redoing fire polygons pmtiles to have each range be its own layer "
|
| 923 |
-
]
|
| 924 |
-
},
|
| 925 |
-
{
|
| 926 |
-
"cell_type": "code",
|
| 927 |
-
"execution_count": null,
|
| 928 |
-
"id": "2161c50b-0328-474f-aa57-215e14fe33c2",
|
| 929 |
-
"metadata": {},
|
| 930 |
-
"outputs": [],
|
| 931 |
-
"source": [
|
| 932 |
-
"def generate_pmtiles(input_file1, input_file2, input_file3, input_file4, output_file, max_zoom=12):\n",
|
| 933 |
-
" # Ensure Tippecanoe is installed\n",
|
| 934 |
-
" if subprocess.call([\"which\", \"tippecanoe\"], stdout=subprocess.DEVNULL) != 0:\n",
|
| 935 |
-
" raise RuntimeError(\"Tippecanoe is not installed or not in PATH\")\n",
|
| 936 |
-
"\n",
|
| 937 |
-
" # Construct the Tippecanoe command\n",
|
| 938 |
-
" command = [\n",
|
| 939 |
-
" \"tippecanoe\",\n",
|
| 940 |
-
" \"-o\", output_file,\n",
|
| 941 |
-
" \"-zg\",\n",
|
| 942 |
-
" \"--extend-zooms-if-still-dropping\",\n",
|
| 943 |
-
" \"--force\",\n",
|
| 944 |
-
" \"--projection\", \"EPSG:4326\", \n",
|
| 945 |
-
" \"-L\",\"layer1:\"+input_file1,\n",
|
| 946 |
-
" \"-L\",\"layer2:\"+input_file2,\n",
|
| 947 |
-
" \"-L\",\"layer3:\"+input_file3,\n",
|
| 948 |
-
" \"-L\",\"layer4:\"+input_file4,\n",
|
| 949 |
-
"\n",
|
| 950 |
-
" ]\n",
|
| 951 |
-
" # Run Tippecanoe\n",
|
| 952 |
-
" try:\n",
|
| 953 |
-
" subprocess.run(command, check=True)\n",
|
| 954 |
-
" print(f\"Successfully generated PMTiles file: {output_file}\")\n",
|
| 955 |
-
" except subprocess.CalledProcessError as e:\n",
|
| 956 |
-
" print(f\"Error running Tippecanoe: {e}\")\n"
|
| 957 |
-
]
|
| 958 |
-
},
|
| 959 |
-
{
|
| 960 |
-
"cell_type": "code",
|
| 961 |
-
"execution_count": null,
|
| 962 |
-
"id": "3a15d11f-ef32-4af3-8b72-b43acd43cf08",
|
| 963 |
-
"metadata": {},
|
| 964 |
-
"outputs": [],
|
| 965 |
-
"source": [
|
| 966 |
-
"rxburn_20 = (con\n",
|
| 967 |
-
" .read_parquet(\"rxburn22_1.parquet\")\n",
|
| 968 |
-
" .rename(year = \"YEAR_\")\n",
|
| 969 |
-
" .filter(_.STATE == \"CA\", _.year != ' ', _.year != '')\n",
|
| 970 |
-
" .cast({\"year\":\"int\"})\n",
|
| 971 |
-
" .filter(_.year>=2003)\n",
|
| 972 |
-
" .mutate(\n",
|
| 973 |
-
" geometry=ibis.ifelse(\n",
|
| 974 |
-
" _.geometry.is_valid(),\n",
|
| 975 |
-
" _.geometry, # Keep the geometry if it's valid\n",
|
| 976 |
-
" _.geometry.buffer(0) # Apply buffer(0) to fix invalid geometries\n",
|
| 977 |
-
" )\n",
|
| 978 |
-
" )\n",
|
| 979 |
-
" .mutate(geometry = _.geometry.convert(\"EPSG:3310\",\"EPSG:4326\"))\n",
|
| 980 |
-
" )\n",
|
| 981 |
-
"\n",
|
| 982 |
-
"rxburn_10 = (rxburn_20.filter(_.year>=2013))\n",
|
| 983 |
-
"rxburn_5 = (rxburn_20.filter(_.year>=2018))\n",
|
| 984 |
-
"rxburn_2 = (rxburn_20.filter(_.year>=2022))\n",
|
| 985 |
-
"\n",
|
| 986 |
-
"rxburn_20_df = rxburn_20.execute().set_crs(\"EPSG:4326\").to_file(\"rxburn_20.geojson\")\n",
|
| 987 |
-
"rxburn_10_df = rxburn_10.execute().set_crs(\"EPSG:4326\").to_file(\"rxburn_10.geojson\")\n",
|
| 988 |
-
"rxburn_5_df = rxburn_5.execute().set_crs(\"EPSG:4326\").to_file(\"rxburn_5.geojson\")\n",
|
| 989 |
-
"rxburn_2_df = rxburn_2.execute().set_crs(\"EPSG:4326\").to_file(\"rxburn_2.geojson\")\n",
|
| 990 |
-
"\n",
|
| 991 |
-
"\n",
|
| 992 |
-
"generate_pmtiles(\"rxburn_20.geojson\",\"rxburn_10.geojson\",\"rxburn_5.geojson\",\"rxburn_2.geojson\",\"cal_rxburn_2022.pmtiles\")\n",
|
| 993 |
-
"hf_upload(\"cal_rxburn_2022.pmtiles\", \"boettiger-lab/ca-30x30\",\"dataset\")\n"
|
| 994 |
-
]
|
| 995 |
-
},
|
| 996 |
-
{
|
| 997 |
-
"cell_type": "code",
|
| 998 |
-
"execution_count": null,
|
| 999 |
-
"id": "1220c348-c68b-4475-ba0f-ef563fea7345",
|
| 1000 |
-
"metadata": {},
|
| 1001 |
-
"outputs": [],
|
| 1002 |
-
"source": [
|
| 1003 |
-
"fire_20 = (con\n",
|
| 1004 |
-
" .read_parquet(\"firep22_1.parquet\")\n",
|
| 1005 |
-
" .rename(year = \"YEAR_\")\n",
|
| 1006 |
-
" .filter(_.STATE == \"CA\", _.year != '')\n",
|
| 1007 |
-
" .cast({\"year\":\"int\"})\n",
|
| 1008 |
-
" .filter(_.year>=2003)\n",
|
| 1009 |
-
" .select(\"year\",\"geometry\")\n",
|
| 1010 |
-
" .mutate(\n",
|
| 1011 |
-
" geometry=ibis.ifelse(\n",
|
| 1012 |
-
" _.geometry.is_valid(),\n",
|
| 1013 |
-
" _.geometry, # Keep the geometry if it's valid\n",
|
| 1014 |
-
" _.geometry.buffer(0) # Apply buffer(0) to fix invalid geometries\n",
|
| 1015 |
-
" )\n",
|
| 1016 |
-
" )\n",
|
| 1017 |
-
" .mutate(geometry = _.geometry.convert(\"EPSG:3310\",\"EPSG:4326\"))\n",
|
| 1018 |
-
" )\n",
|
| 1019 |
-
"\n",
|
| 1020 |
-
"fire_10 = (fire_20.filter(_.year>=2013))\n",
|
| 1021 |
-
"fire_5 = (fire_20.filter(_.year>=2018))\n",
|
| 1022 |
-
"fire_2 = (fire_20.filter(_.year>=2022))\n",
|
| 1023 |
-
"\n",
|
| 1024 |
-
"fire_20_df = fire_20.execute().set_crs(\"EPSG:4326\").to_file(\"fire_20.geojson\")\n",
|
| 1025 |
-
"fire_10_df = fire_10.execute().set_crs(\"EPSG:4326\").to_file(\"fire_10.geojson\")\n",
|
| 1026 |
-
"fire_5_df = fire_5.execute().set_crs(\"EPSG:4326\").to_file(\"fire_5.geojson\")\n",
|
| 1027 |
-
"fire_2_df = fire_2.execute().set_crs(\"EPSG:4326\").to_file(\"fire_2.geojson\")\n",
|
| 1028 |
-
"\n",
|
| 1029 |
-
"\n",
|
| 1030 |
-
"generate_pmtiles(\"fire_20.geojson\",\"fire_10.geojson\",\"fire_5.geojson\",\"fire_2.geojson\",\"cal_fire_2022.pmtiles\")\n",
|
| 1031 |
-
"hf_upload(\"cal_fire_2022.pmtiles\", \"boettiger-lab/ca-30x30\",\"dataset\")\n"
|
| 1032 |
-
]
|
| 1033 |
-
}
|
| 1034 |
-
],
|
| 1035 |
-
"metadata": {
|
| 1036 |
-
"kernelspec": {
|
| 1037 |
-
"display_name": "Python 3 (ipykernel)",
|
| 1038 |
-
"language": "python",
|
| 1039 |
-
"name": "python3"
|
| 1040 |
-
},
|
| 1041 |
-
"language_info": {
|
| 1042 |
-
"codemirror_mode": {
|
| 1043 |
-
"name": "ipython",
|
| 1044 |
-
"version": 3
|
| 1045 |
-
},
|
| 1046 |
-
"file_extension": ".py",
|
| 1047 |
-
"mimetype": "text/x-python",
|
| 1048 |
-
"name": "python",
|
| 1049 |
-
"nbconvert_exporter": "python",
|
| 1050 |
-
"pygments_lexer": "ipython3",
|
| 1051 |
-
"version": "3.12.7"
|
| 1052 |
-
}
|
| 1053 |
-
},
|
| 1054 |
-
"nbformat": 4,
|
| 1055 |
-
"nbformat_minor": 5
|
| 1056 |
-
}
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|