{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "\n", "import pystac_client\n", "import planetary_computer\n", "import odc.stac\n", "import geopandas as gpd\n", "import dask.distributed\n", "import matplotlib.pyplot as plt\n", "import rioxarray\n", "from datetime import datetime, timedelta\n", "\n", "\n", "\n", "nps = gpd.read_file(\"/vsicurl/https://huggingface.co/datasets/cboettig/biodiversity/resolve/main/data/NPS.gdb\")\n", "calfire = gpd.read_file(\"/vsicurl/https://huggingface.co/datasets/cboettig/biodiversity/resolve/main/data/fire22_1.gdb\", layer = \"firep22_1\")\n", "# fire = gpd.read_file(\"/vsizip/vsicurl/https://edcintl.cr.usgs.gov/downloads/sciweb1/shared/MTBS_Fire/data/composite_data/burned_area_extent_shapefile/mtbs_perimeter_data.zip\"\n", "\n", "# extract and reproject the Joshua Tree NP Polygon\n", "jtree = nps[nps.PARKNAME == \"Joshua Tree\"].to_crs(calfire.crs)\n", "\n", "# All Fires in the DB that intersect the Park\n", "jtree_fires = jtree.overlay(calfire, how=\"intersection\")\n", "\n", "# Extract a polygon if interest. > 2015 for Sentinel, otherwise we can use LandSat\n", "recent = jtree_fires[jtree_fires.YEAR_ > \"2015\"]\n", "big = recent[recent.Shape_Area == recent.Shape_Area.max()].to_crs(\"EPSG:4326\")\n", "\n", "# Get bounding box + dates before & after fire for STAC search\n", "box = big.buffer(0.01).bounds.to_numpy()[0] # Fire bbox + buffer\n", "alarm_date = datetime.strptime(big.ALARM_DATE.item(), \"%Y-%m-%dT%H:%M:%S+00:00\") \n", "before_date = alarm_date - timedelta(days=14)\n", "after_date = alarm_date + timedelta(days=14)\n", "search_dates = before_date.strftime(\"%Y-%m-%d\") + \"/\" + after_date.strftime(\"%Y-%m-%d\")\n", "\n" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "\n", "def stac_search(box, datetime): \n", " # STAC Search for this imagery in space/time window\n", " items = (\n", " pystac_client.Client.\n", " open(\"https://cmr.earthdata.nasa.gov/stac/LPCLOUD\").\n", " search(collections=[\"HLSL30.v2.0\"],\n", " bbox=box,\n", " datetime=datetime).\n", " item_collection())\n", " return items\n", "\n", "def compute_ndvi(items, box):\n", " # Time to compute:\n", " client = dask.distributed.Client()\n", " # landsat_bands = [\"nir08\", \"swir16\"]\n", " sentinel_bands = [\"B05\", \"B04\", \"Fmask\"] # NIR, SWIR, and Cloud Mask\n", "\n", " # The magic of gdalwarper. Can also resample, reproject, and aggregate on the fly\n", " data = odc.stac.load(items,\n", " bands=sentinel_bands,\n", " bbox=box,\n", " crs=\"EPSG:32611\"\n", " )\n", "\n", " red = data[\"B04\"].astype(\"float\")\n", " nir = data[\"B05\"].astype(\"float\")\n", "\n", " # can resample and aggregate in xarray. compute with dask\n", " ndvi = ((nir - red) / (nir + red)).compute()\n", " \n", " return ndvi\n" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "
<xarray.DataArray 'B04' (time: 7, latitude: 1, longitude: 1)>\n",
"array([[[nan]],\n",
"\n",
" [[nan]],\n",
"\n",
" [[nan]],\n",
"\n",
" [[nan]],\n",
"\n",
" [[nan]],\n",
"\n",
" [[nan]],\n",
"\n",
" [[nan]]], dtype=float32)\n",
"Coordinates:\n",
" * latitude (latitude) float64 45.0\n",
" * longitude (longitude) float64 -105.0\n",
" spatial_ref int32 4326\n",
" * time (time) datetime64[ns] 2022-05-15T18:22:05.159000 ... 2022-06...