|
| 1 | +{ |
| 2 | + "cells": [ |
| 3 | + { |
| 4 | + "cell_type": "markdown", |
| 5 | + "id": "ccb0fe8e-6b86-403f-84e0-c910f70bfc23", |
| 6 | + "metadata": {}, |
| 7 | + "source": [ |
| 8 | + "# JupyterGIS demo" |
| 9 | + ] |
| 10 | + }, |
| 11 | + { |
| 12 | + "cell_type": "markdown", |
| 13 | + "id": "6a1dcd59-5bd5-4e1a-93d3-8e0850aacc9e", |
| 14 | + "metadata": {}, |
| 15 | + "source": [ |
| 16 | + "## Outline\n", |
| 17 | + "\n", |
| 18 | + "* Aggregate gridded data based on vector regions (e.g. neighborhoods) \n", |
| 19 | + " * Not straightforward to do in Python \n", |
| 20 | + " * Design:\n", |
| 21 | + " * Start in a Notebook, prepared with Maryam’s expertise \n", |
| 22 | + " * Loading GeoPandas, tools for Zonal Statistics \n", |
| 23 | + " * Programmatically create .jGIS document, add input data sources and output data sources. \n", |
| 24 | + " * Demonstrate collaboration of JGIS alongside Notebook. Annotation, ad layer from catalog, etc. \n", |
| 25 | + " * Examples: \n", |
| 26 | + " * [https://medium.com/data-science/zonal-statistics-algorithm-with-python-in-4-steps-382a3b66648a](https://medium.com/data-science/zonal-statistics-algorithm-with-python-in-4-steps-382a3b66648a) \n", |
| 27 | + " * [https://automating-gis-processes.github.io/CSC18/lessons/L6/zonal-statistics.html](https://automating-gis-processes.github.io/CSC18/lessons/L6/zonal-statistics.html) \n", |
| 28 | + " * [https://carpentries-incubator.github.io/geospatial-python/10-zonal-statistics.html](https://carpentries-incubator.github.io/geospatial-python/10-zonal-statistics.html) \n", |
| 29 | + " * From Carl: [https://espm-288.carlboettiger.info/tutorials/python/spatial-4.html](https://espm-288.carlboettiger.info/tutorials/python/spatial-4.html)" |
| 30 | + ] |
| 31 | + }, |
| 32 | + { |
| 33 | + "cell_type": "markdown", |
| 34 | + "id": "0b259e2b-0758-430d-ba4e-d48184395473", |
| 35 | + "metadata": {}, |
| 36 | + "source": [ |
| 37 | + "## From geopythontutorials.com\n", |
| 38 | + "\n", |
| 39 | + "https://www.geopythontutorials.com/notebooks/xarray_zonal_stats.html?utm_source=chatgpt.com\n", |
| 40 | + "\n", |
| 41 | + "New dependencies\n", |
| 42 | + "\n", |
| 43 | + "* rioxarray\n", |
| 44 | + "* geocube\n", |
| 45 | + "* xarray-spatial" |
| 46 | + ] |
| 47 | + }, |
| 48 | + { |
| 49 | + "cell_type": "markdown", |
| 50 | + "id": "a2a3aa9d-1a0e-4065-b7c3-b99492a81800", |
| 51 | + "metadata": {}, |
| 52 | + "source": [ |
| 53 | + "### Download the data" |
| 54 | + ] |
| 55 | + }, |
| 56 | + { |
| 57 | + "cell_type": "code", |
| 58 | + "execution_count": null, |
| 59 | + "id": "6b7b45a4-c41a-4c68-8fa2-d7472898af0a", |
| 60 | + "metadata": {}, |
| 61 | + "outputs": [], |
| 62 | + "source": [ |
| 63 | + "import os\n", |
| 64 | + "\n", |
| 65 | + "data_folder = \"data\"\n", |
| 66 | + "\n", |
| 67 | + "def download(url, data_folder):\n", |
| 68 | + " filename = os.path.join(data_folder, os.path.basename(url))\n", |
| 69 | + " if not os.path.exists(filename):\n", |
| 70 | + " from urllib.request import urlretrieve\n", |
| 71 | + " local, _ = urlretrieve(url, filename)\n", |
| 72 | + " print('Downloaded ' + local)\n", |
| 73 | + "\n", |
| 74 | + "raster_file = 'chirps-v2.0.2021.tif'\n", |
| 75 | + "zones_file = 'cb_2021_us_county_500k.zip'\n", |
| 76 | + "\n", |
| 77 | + "files = [\n", |
| 78 | + " 'https://data.chc.ucsb.edu/products/CHIRPS-2.0/global_annual/tifs/' + raster_file,\n", |
| 79 | + " 'https://www2.census.gov/geo/tiger/GENZ2021/shp/' + zones_file,\n", |
| 80 | + "]\n", |
| 81 | + "\n", |
| 82 | + "for file in files:\n", |
| 83 | + " download(file, data_folder)" |
| 84 | + ] |
| 85 | + }, |
| 86 | + { |
| 87 | + "cell_type": "markdown", |
| 88 | + "id": "742d07d7-e745-4e9f-8c46-2662be64e54f", |
| 89 | + "metadata": {}, |
| 90 | + "source": [ |
| 91 | + "### Data pre-processing" |
| 92 | + ] |
| 93 | + }, |
| 94 | + { |
| 95 | + "cell_type": "code", |
| 96 | + "execution_count": null, |
| 97 | + "id": "17aa7d8f-51a5-4022-b784-5c6021053f65", |
| 98 | + "metadata": {}, |
| 99 | + "outputs": [], |
| 100 | + "source": [ |
| 101 | + "import geopandas as gpd\n", |
| 102 | + "\n", |
| 103 | + "zones_file_path = os.path.join(data_folder, zones_file)\n", |
| 104 | + "\n", |
| 105 | + "zones_df = gpd.read_file(zones_file_path)\n", |
| 106 | + "# TODO: Louisiana instead?\n", |
| 107 | + "california_df = zones_df[zones_df['STATE_NAME'] == 'California'].copy()\n", |
| 108 | + "california_df.iloc[:5, :5]" |
| 109 | + ] |
| 110 | + }, |
| 111 | + { |
| 112 | + "cell_type": "code", |
| 113 | + "execution_count": null, |
| 114 | + "id": "68ead58b-fd90-4a30-8cc3-f94b362f09f0", |
| 115 | + "metadata": {}, |
| 116 | + "outputs": [], |
| 117 | + "source": [ |
| 118 | + "california_df['GEOID'] = california_df.GEOID.astype(int)" |
| 119 | + ] |
| 120 | + }, |
| 121 | + { |
| 122 | + "cell_type": "code", |
| 123 | + "execution_count": null, |
| 124 | + "id": "da98d0a2-b1cc-4a6b-8754-2df48918c3b9", |
| 125 | + "metadata": {}, |
| 126 | + "outputs": [], |
| 127 | + "source": [ |
| 128 | + "import rioxarray as rxr\n", |
| 129 | + "\n", |
| 130 | + "raster_filepath = os.path.join(data_folder, raster_file)\n", |
| 131 | + "raster = rxr.open_rasterio(raster_filepath, mask_and_scale=True)\n", |
| 132 | + "clipped = raster.rio.clip(california_df.geometry)\n", |
| 133 | + "clipped" |
| 134 | + ] |
| 135 | + }, |
| 136 | + { |
| 137 | + "cell_type": "code", |
| 138 | + "execution_count": null, |
| 139 | + "id": "e2f81cff-8836-43c1-9477-4950b771bc79", |
| 140 | + "metadata": {}, |
| 141 | + "outputs": [], |
| 142 | + "source": [ |
| 143 | + "precipitation = clipped.sel(band=1)\n", |
| 144 | + "precipitation" |
| 145 | + ] |
| 146 | + }, |
| 147 | + { |
| 148 | + "cell_type": "code", |
| 149 | + "execution_count": null, |
| 150 | + "id": "0b1f55b7-081b-43ec-9bbb-55295a1cb42a", |
| 151 | + "metadata": {}, |
| 152 | + "outputs": [], |
| 153 | + "source": [ |
| 154 | + "from geocube.api.core import make_geocube\n", |
| 155 | + "\n", |
| 156 | + "california_raster = make_geocube(\n", |
| 157 | + " vector_data=california_df,\n", |
| 158 | + " measurements=['GEOID'],\n", |
| 159 | + " like=precipitation,\n", |
| 160 | + ")\n", |
| 161 | + "california_raster" |
| 162 | + ] |
| 163 | + }, |
| 164 | + { |
| 165 | + "cell_type": "code", |
| 166 | + "execution_count": null, |
| 167 | + "id": "216b6a74-b543-4896-8348-f84b80babecb", |
| 168 | + "metadata": {}, |
| 169 | + "outputs": [], |
| 170 | + "source": [ |
| 171 | + "from xrspatial import zonal_stats\n", |
| 172 | + "\n", |
| 173 | + "stats_df = zonal_stats(zones=california_raster.GEOID, values=precipitation)\n", |
| 174 | + "stats_df.iloc[:5]" |
| 175 | + ] |
| 176 | + }, |
| 177 | + { |
| 178 | + "cell_type": "code", |
| 179 | + "execution_count": null, |
| 180 | + "id": "ae05145f-b31c-4d10-8c3e-852c722fa0ec", |
| 181 | + "metadata": {}, |
| 182 | + "outputs": [], |
| 183 | + "source": [ |
| 184 | + "stats_df['GEOID'] = stats_df['zone'].astype(int)" |
| 185 | + ] |
| 186 | + }, |
| 187 | + { |
| 188 | + "cell_type": "code", |
| 189 | + "execution_count": null, |
| 190 | + "id": "aa371396-e42b-4c18-a07f-7d2eb8cb0ad6", |
| 191 | + "metadata": {}, |
| 192 | + "outputs": [], |
| 193 | + "source": [ |
| 194 | + "joined = california_df.merge(stats_df[['GEOID', 'mean']], on='GEOID')\n", |
| 195 | + "joined.iloc[:5, -5:]" |
| 196 | + ] |
| 197 | + }, |
| 198 | + { |
| 199 | + "cell_type": "code", |
| 200 | + "execution_count": null, |
| 201 | + "id": "e306fe32-9197-455a-aedb-0d3b24fb215d", |
| 202 | + "metadata": {}, |
| 203 | + "outputs": [], |
| 204 | + "source": [ |
| 205 | + "import matplotlib.pyplot as plt\n", |
| 206 | + "\n", |
| 207 | + "fig, ax = plt.subplots(1, 1)\n", |
| 208 | + "fig.set_size_inches(10,10)\n", |
| 209 | + "\n", |
| 210 | + "legend_kwds={\n", |
| 211 | + " 'orientation': 'horizontal', # Make the legend horizontal\n", |
| 212 | + " 'shrink': 0.5, # Reduce the size of the legend bar by 50%\n", |
| 213 | + " 'pad': 0.05, # Add some padding around the legend\n", |
| 214 | + " 'label': 'Precipitation (mm)', # Set the legend label (optional)\n", |
| 215 | + " }\n", |
| 216 | + "joined.plot(ax=ax, column='mean', cmap='Blues',\n", |
| 217 | + " legend=True, legend_kwds=legend_kwds)\n", |
| 218 | + "ax.set_axis_off()\n", |
| 219 | + "ax.set_title('Total Precipitation 2021 for California Counties')\n", |
| 220 | + "plt.show()" |
| 221 | + ] |
| 222 | + }, |
| 223 | + { |
| 224 | + "cell_type": "code", |
| 225 | + "execution_count": null, |
| 226 | + "id": "0b1d1499-fb42-437e-a979-88ddb74f4822", |
| 227 | + "metadata": {}, |
| 228 | + "outputs": [], |
| 229 | + "source": [ |
| 230 | + "joined.explore()" |
| 231 | + ] |
| 232 | + }, |
| 233 | + { |
| 234 | + "cell_type": "markdown", |
| 235 | + "id": "8d63cc59-78cd-4639-9c96-f6ac3e0dd208", |
| 236 | + "metadata": { |
| 237 | + "jp-MarkdownHeadingCollapsed": true |
| 238 | + }, |
| 239 | + "source": [ |
| 240 | + "## From Carl's class" |
| 241 | + ] |
| 242 | + }, |
| 243 | + { |
| 244 | + "cell_type": "markdown", |
| 245 | + "id": "64b27ffc-7f33-4374-b0a8-160bd3c50567", |
| 246 | + "metadata": {}, |
| 247 | + "source": [ |
| 248 | + "New dependencies (don't add to environment, this is just for accessing data):\n", |
| 249 | + "\n", |
| 250 | + "* ibis-duckdb\n", |
| 251 | + "* odc\n", |
| 252 | + "\n", |
| 253 | + "New dependencies (add to environment):\n", |
| 254 | + "\n", |
| 255 | + "* exactextract" |
| 256 | + ] |
| 257 | + }, |
| 258 | + { |
| 259 | + "cell_type": "markdown", |
| 260 | + "id": "936e26c2-05ca-48b4-990d-7097b51e93ac", |
| 261 | + "metadata": {}, |
| 262 | + "source": [ |
| 263 | + "### Setting up raster data (NDVI)\n", |
| 264 | + "\n", |
| 265 | + "https://espm-288.carlboettiger.info/tutorials/python/spatial-3.html" |
| 266 | + ] |
| 267 | + }, |
| 268 | + { |
| 269 | + "cell_type": "code", |
| 270 | + "execution_count": null, |
| 271 | + "id": "21dc77b3-a2dc-4356-86ff-6a9d617231d2", |
| 272 | + "metadata": { |
| 273 | + "jupyter": { |
| 274 | + "source_hidden": true |
| 275 | + } |
| 276 | + }, |
| 277 | + "outputs": [], |
| 278 | + "source": [ |
| 279 | + "!export CURL_CA_BUNDLE=/etc/ssl/certs/ca-certificates.crt" |
| 280 | + ] |
| 281 | + }, |
| 282 | + { |
| 283 | + "cell_type": "code", |
| 284 | + "execution_count": null, |
| 285 | + "id": "edf53f2b-657d-4683-9518-8257a93e5a21", |
| 286 | + "metadata": { |
| 287 | + "jupyter": { |
| 288 | + "source_hidden": true |
| 289 | + } |
| 290 | + }, |
| 291 | + "outputs": [], |
| 292 | + "source": [ |
| 293 | + "import ibis\n", |
| 294 | + "from ibis import _\n", |
| 295 | + "\n", |
| 296 | + "con = ibis.duckdb.connect(extensions=[\"spatial\"])\n", |
| 297 | + "# TODO: Why doesn't this work?\n", |
| 298 | + "# \"does not exist in the file system, and is not recognized as a supported dataset name\"\n", |
| 299 | + "#redlines = (\n", |
| 300 | + "# con\n", |
| 301 | + "# .read_geo(\"/vsicurl/https://dsl.richmond.edu/panorama/redlining/static/mappinginequality.gpkg\")\n", |
| 302 | + "# .filter(_.city == \"New Haven\", _.residential)\n", |
| 303 | + "#)\n", |
| 304 | + "redlines = (\n", |
| 305 | + " con\n", |
| 306 | + " .read_geo(\"./mappinginequality.gpkg\")\n", |
| 307 | + " .filter(_.city == \"New Haven\", _.residential)\n", |
| 308 | + ")\n", |
| 309 | + "city = redlines.execute()\n", |
| 310 | + "box = city.total_bounds\n", |
| 311 | + "box" |
| 312 | + ] |
| 313 | + }, |
| 314 | + { |
| 315 | + "cell_type": "code", |
| 316 | + "execution_count": null, |
| 317 | + "id": "61363fcb-e942-489e-a631-9cfb5db72800", |
| 318 | + "metadata": { |
| 319 | + "jupyter": { |
| 320 | + "source_hidden": true |
| 321 | + } |
| 322 | + }, |
| 323 | + "outputs": [], |
| 324 | + "source": [ |
| 325 | + "from pystac_client import Client\n", |
| 326 | + "\n", |
| 327 | + "items = (\n", |
| 328 | + " Client.\n", |
| 329 | + " open(\"https://earth-search.aws.element84.com/v1\").\n", |
| 330 | + " search(\n", |
| 331 | + " collections = ['sentinel-2-l2a'],\n", |
| 332 | + " bbox = box,\n", |
| 333 | + " datetime = \"2024-06-01/2024-09-01\",\n", |
| 334 | + " query={\"eo:cloud_cover\": {\"lt\": 20}}).\n", |
| 335 | + " item_collection()\n", |
| 336 | + ")\n", |
| 337 | + "items" |
| 338 | + ] |
| 339 | + }, |
| 340 | + { |
| 341 | + "cell_type": "code", |
| 342 | + "execution_count": null, |
| 343 | + "id": "814ce3ab-c48b-4ee5-b7b1-2572ee6fb305", |
| 344 | + "metadata": { |
| 345 | + "jupyter": { |
| 346 | + "source_hidden": true |
| 347 | + } |
| 348 | + }, |
| 349 | + "outputs": [], |
| 350 | + "source": [ |
| 351 | + "import odc.stac" |
| 352 | + ] |
| 353 | + }, |
| 354 | + { |
| 355 | + "cell_type": "code", |
| 356 | + "execution_count": null, |
| 357 | + "id": "50adaad1-493a-47e2-843c-fa596000dfac", |
| 358 | + "metadata": {}, |
| 359 | + "outputs": [], |
| 360 | + "source": [ |
| 361 | + "### Zonal statistics\n", |
| 362 | + "\n", |
| 363 | + "https://espm-288.carlboettiger.info/tutorials/python/spatial-4.html" |
| 364 | + ] |
| 365 | + } |
| 366 | + ], |
| 367 | + "metadata": { |
| 368 | + "kernelspec": { |
| 369 | + "display_name": "Python 3 (ipykernel)", |
| 370 | + "language": "python", |
| 371 | + "name": "python3" |
| 372 | + }, |
| 373 | + "language_info": { |
| 374 | + "codemirror_mode": { |
| 375 | + "name": "ipython", |
| 376 | + "version": 3 |
| 377 | + }, |
| 378 | + "file_extension": ".py", |
| 379 | + "mimetype": "text/x-python", |
| 380 | + "name": "python", |
| 381 | + "nbconvert_exporter": "python", |
| 382 | + "pygments_lexer": "ipython3", |
| 383 | + "version": "3.13.0" |
| 384 | + } |
| 385 | + }, |
| 386 | + "nbformat": 4, |
| 387 | + "nbformat_minor": 5 |
| 388 | +} |
0 commit comments