|
1 | 1 | { |
2 | 2 | "cells": [ |
3 | 3 | { |
4 | | - "attachments": {}, |
5 | 4 | "cell_type": "markdown", |
6 | 5 | "metadata": {}, |
7 | 6 | "source": [ |
|
10 | 9 | ] |
11 | 10 | }, |
12 | 11 | { |
13 | | - "attachments": {}, |
14 | 12 | "cell_type": "markdown", |
15 | 13 | "metadata": {}, |
16 | 14 | "source": [ |
17 | 15 | "## Overview\n", |
18 | 16 | "\n", |
19 | | - "In this notebook, we will take a look at how to retrieve Sentinel-2 L2A satellite imagery from the [Microsoft Planetary Computer Data Catalog (MSPC)](https://planetarycomputer.microsoft.com/catalog). We will go over how to interact with the Data Catalog, which exposes a [SpatioTemporal Asset Catalog (STAC)](https://stacspec.org/en) interface for querying, searching and retrieving data. We will use the [stackstac](https://stackstac.readthedocs.io/en/latest/) package to load the data lazily, which means data is not *actually* read unless required (say, for plotting). Once loaded, we will process the data and make a simple interactive dashboard to look at the satellite imagery over a location for different seasons. We will use the [HoloViz ecosystem](https://holoviz.org/background.html) for the interactive dashboard." |
| 17 | + "In this notebook, we will take a look at how to retrieve Sentinel-2 L2A satellite imagery from the [Microsoft Planetary Computer Data Catalog (MSPC)](https://planetarycomputer.microsoft.com/catalog). We will go over how to interact with the Data Catalog, which exposes a [SpatioTemporal Asset Catalog (STAC)](https://stacspec.org/en) interface for querying, searching and retrieving data. We will use the [stackstac](https://stackstac.readthedocs.io/en/latest/) package to load the data lazily, which means data is not *actually* read unless required (say, for plotting). Once loaded, we will process the data and make a simple interactive dashboard to look at the satellite imagery over a location for different seasons. We will use the [HoloViz ecosystem](https://holoviz.org/background.html) for the interactive dashboard.\n", |
| 18 | + "\n", |
| 19 | + "# TODO: Add authorship using CITATION.cff" |
20 | 20 | ] |
21 | 21 | }, |
22 | 22 | { |
|
48 | 48 | { |
49 | 49 | "cell_type": "code", |
50 | 50 | "execution_count": null, |
51 | | - "metadata": {}, |
| 51 | + "metadata": { |
| 52 | + "tags": [] |
| 53 | + }, |
52 | 54 | "outputs": [], |
53 | 55 | "source": [ |
54 | 56 | "import os\n", |
55 | 57 | "import pandas as pd\n", |
| 58 | + "import numpy as np\n", |
56 | 59 | "import xarray as xr\n", |
57 | 60 | "import stackstac\n", |
58 | 61 | "import pystac_client\n", |
|
65 | 68 | "from pystac.extensions.eo import EOExtension as eo\n", |
66 | 69 | "import datetime\n", |
67 | 70 | "from cartopy import crs\n", |
| 71 | + "import dask\n", |
68 | 72 | "from dask.distributed import Client, LocalCluster\n", |
| 73 | + "import odc.stac\n", |
69 | 74 | "\n", |
70 | 75 | "xr.set_options(keep_attrs=True)\n", |
71 | 76 | "hv.extension('bokeh')\n", |
|
91 | 96 | ] |
92 | 97 | }, |
93 | 98 | { |
94 | | - "attachments": {}, |
95 | 99 | "cell_type": "markdown", |
96 | 100 | "metadata": {}, |
97 | 101 | "source": [ |
|
144 | 148 | ] |
145 | 149 | }, |
146 | 150 | { |
147 | | - "attachments": {}, |
148 | 151 | "cell_type": "markdown", |
149 | 152 | "metadata": {}, |
150 | 153 | "source": [ |
|
173 | 176 | "bbox = [-105.283263,39.972809,-105.266569,39.987640] # NCAR, boulder, CO. bbox from http://bboxfinder.com/\n", |
174 | 177 | "date_range = \"2022-01-01/2022-12-31\"\n", |
175 | 178 | "collection = \"sentinel-2-l2a\" # full id of collection\n", |
176 | | - "cloud_thresh = 40" |
| 179 | + "cloud_thresh = 30" |
177 | 180 | ] |
178 | 181 | }, |
179 | 182 | { |
|
193 | 196 | ] |
194 | 197 | }, |
195 | 198 | { |
196 | | - "attachments": {}, |
197 | 199 | "cell_type": "markdown", |
198 | 200 | "metadata": {}, |
199 | 201 | "source": [ |
|
223 | 225 | "|B03|Green|\n", |
224 | 226 | "|B02|Blue|\n", |
225 | 227 | "\n", |
226 | | - "We will use the `stackstac.stack` function to load in the assets that start with the alphabet 'B'. This function will return a lazily-loaded `xr.DataArray` (using dask). " |
| 228 | + "We will use the `odc.stac.stac_load` function to load in the assets that start with the alphabet 'B'. This function will return a lazily-loaded `xr.DataSet` (using dask). For plotting purposes it is better if we have the data as a `xr.DataArray` instead with the bands as a dimension. We can do that using `.to_array(dim=<dim_name>)` method of a dataset." |
227 | 229 | ] |
228 | 230 | }, |
229 | 231 | { |
|
234 | 236 | "source": [ |
235 | 237 | "bands_of_interest = [b for b in all_bands if b.startswith('B')]\n", |
236 | 238 | "\n", |
237 | | - "da = stackstac.stack(\n", |
| 239 | + "da = odc.stac.stac_load(\n", |
238 | 240 | " items,\n", |
239 | | - " bounds_latlon=bbox,\n", |
240 | | - " assets=bands_of_interest,\n", |
241 | | - " chunksize='50MiB'\n", |
242 | | - ")\n", |
| 241 | + " bands=bands_of_interest,\n", |
| 242 | + " bbox=bbox,\n", |
| 243 | + " chunks={}, # <-- use Dask\n", |
| 244 | + ").to_array(dim='band')\n", |
243 | 245 | "da" |
244 | 246 | ] |
245 | 247 | }, |
|
303 | 305 | "metadata": {}, |
304 | 306 | "source": [ |
305 | 307 | "Now that we have a harmonized dataset, we still need to process the data as follows:\n", |
306 | | - "- On closer inspection, I found some duplicate data in the `time` dimension which were leading to errors later in the notebook. We can drop the duplicate values using the `da.drop_duplicates` method.\n", |
307 | | - "- Sentinel-2 L2A provides the Surface Reflectance (SR) data, which usually ranges from 0 (no reflection) to 1.0 (complete reflection). However, the actual values in the loaded dataset ranges from 0 to ~10,000. These data values need to be scaled to 0.0-1.0 by dividing the data by 10,000. More details can be found in [section 2.3.10 of this document](https://sentinel.esa.int/documents/247904/685211/Sen2-Cor-L2A-Input-Output-Data-Definition-Document.pdf/e2dd6f01-c9c7-494d-a7f2-cd3be9ad891a?t=1506524754000)." |
| 308 | + "- Sentinel-2 L2A provides the Surface Reflectance (SR) data, which usually ranges from 0 (no reflection) to 1.0 (complete reflection). However, the actual values in the loaded dataset ranges from 0 to ~10,000. These data values need to be scaled to 0.0-1.0 by dividing the data by 10,000. More details can be found in [section 2.3.10 of this document](https://sentinel.esa.int/documents/247904/685211/Sen2-Cor-L2A-Input-Output-Data-Definition-Document.pdf/e2dd6f01-c9c7-494d-a7f2-cd3be9ad891a?t=1506524754000).\n", |
| 309 | + "\n", |
| 310 | + "We will then explicitly trigger the dask computation using the `compute()` method and load the result into memory. This is to reduce repeated calls to retrieve data from MSPC. By loading the processed This wouldn't have been possible if the dataset was large." |
308 | 311 | ] |
309 | 312 | }, |
310 | 313 | { |
|
313 | 316 | "metadata": {}, |
314 | 317 | "outputs": [], |
315 | 318 | "source": [ |
316 | | - "# Seems like there is duplicate data in the time dimension\n", |
317 | | - "da = da.drop_duplicates(dim='time')\n", |
318 | 319 | "da = da / 1e4 # Scale data values from 0:10000 to 0:1.0\n", |
319 | | - "da = da / da.max(dim='band') # Stretch to min-max for *crispier* plots\n", |
320 | | - "da" |
| 320 | + "da = da / da.max(dim='band') # additionally scale from 0-max -> 0-1 for visual quality\n", |
| 321 | + "da = da.compute()" |
321 | 322 | ] |
322 | 323 | }, |
323 | 324 | { |
324 | 325 | "cell_type": "markdown", |
325 | 326 | "metadata": {}, |
326 | 327 | "source": [ |
327 | | - "We have now processed the data so that we can visualize it! *Note: The computation has not been done yet, it will be triggered as soon as we plot the data. This is possible because until now, `dask` has only created the \"task graph\" and we have not yet performed any operation that would trigger the computation yet*.\n", |
| 328 | + "We have now processed the data so that we can visualize it!\n", |
328 | 329 | "\n", |
329 | | - "Let's create a function that will take a `time` input and have do the following tasks:\n", |
330 | | - " 1. plot an interactive RGB image of the data and overlay it on a map of the world.\n", |
331 | | - " 2. provide a [date slider widget](https://panel.holoviz.org/reference/widgets/DateSlider.html) which can be used to interact with the plot.\n", |
332 | | - " 3. only set the default value of the date slider to the `time`, but allow the user to slide through the length of the entire dataset." |
| 330 | + "Let's look at the Blue, Green and Red bands." |
333 | 331 | ] |
334 | 332 | }, |
335 | 333 | { |
|
338 | 336 | "metadata": {}, |
339 | 337 | "outputs": [], |
340 | 338 | "source": [ |
341 | | - "season_names = {\n", |
342 | | - " 1: 'Winter',\n", |
343 | | - " 2: 'Spring',\n", |
344 | | - " 3: 'Summer',\n", |
345 | | - " 4: 'Fall'\n", |
346 | | - "}\n", |
| 339 | + "da.sel(band='B04').isel(time=0).hvplot(x='x', y='y', data_aspect=1, cmap='Blues') \\\n", |
| 340 | + "+ da.sel(band='B03').isel(time=0).hvplot(x='x', y='y', data_aspect=1, cmap='Greens') \\\n", |
| 341 | + "+ da.sel(band='B02').isel(time=0).hvplot(x='x', y='y', data_aspect=1, cmap='Reds')" |
| 342 | + ] |
| 343 | + }, |
| 344 | + { |
| 345 | + "cell_type": "markdown", |
| 346 | + "metadata": {}, |
| 347 | + "source": [ |
| 348 | + "Let us make a dashboard composed of 4 different interactive plots showing the RGB view of the satellite observation for four different seasons.\n", |
| 349 | + "We need a function that will take a `time` input and does the following tasks:\n", |
| 350 | + " 1. plot an interactive RGB image of the data and overlay it on a map of the world.\n", |
| 351 | + " 2. provide a [date slider widget](https://panel.holoviz.org/reference/widgets/DateSlider.html) which can be used to interact with the plot.\n", |
| 352 | + " 3. only set the default value of the date slider to the `time`, but allow the user to slide through the length of the entire dataset.\n", |
347 | 353 | "\n", |
| 354 | + "Using this function, we will be able to compose the dashboard." |
| 355 | + ] |
| 356 | + }, |
| 357 | + { |
| 358 | + "cell_type": "code", |
| 359 | + "execution_count": null, |
| 360 | + "metadata": { |
| 361 | + "tags": [] |
| 362 | + }, |
| 363 | + "outputs": [], |
| 364 | + "source": [ |
348 | 365 | "def rgb_during(time):\n", |
| 366 | + " season_names = {\n", |
| 367 | + " 1: 'Winter',\n", |
| 368 | + " 2: 'Spring',\n", |
| 369 | + " 3: 'Summer',\n", |
| 370 | + " 4: 'Fall'\n", |
| 371 | + " }\n", |
349 | 372 | " da_rgb = da.sel(band=['B04', 'B03', 'B02'])\n", |
350 | 373 | " start_date = pd.to_datetime(da_rgb['time'].min().data).to_pydatetime()\n", |
351 | 374 | " end_date = pd.to_datetime(da_rgb['time'].max().data).to_pydatetime()\n", |
|
355 | 378 | " def get_obs_on(t):\n", |
356 | 379 | " season_key = [month%12 // 3 + 1 for month in range(1, 13)][t.month-1]\n", |
357 | 380 | " season = season_names[season_key]\n", |
358 | | - " return da_rgb.sel(time=t, method='nearest').hvplot.rgb(x='x', y='y', bands='band', data_aspect=1, geo=True, tiles='ESRI', rasterize=True, title=f\"{season}: {t.strftime('%Y-%m-%d')}\")\n", |
| 381 | + " return da.sel(band=['B04', 'B03', 'B02']).sel(time=t, method='nearest').transpose('y', 'x', 'band').hvplot.rgb(\n", |
| 382 | + " x='x', y='y', bands='band', \n", |
| 383 | + " geo=True, tiles='ESRI', crs=crs.epsg(items[0].properties['proj:epsg']), \n", |
| 384 | + " rasterize=True, title=f\"{season}: {t.strftime('%Y-%m-%d')}\")\n", |
| 385 | + " \n", |
359 | 386 | " \n", |
360 | 387 | " return pn.panel(pn.Column(\n", |
361 | 388 | " pn.bind(get_obs_on, t=dt_slider), \n", |
|
366 | 393 | " ))" |
367 | 394 | ] |
368 | 395 | }, |
| 396 | + { |
| 397 | + "cell_type": "code", |
| 398 | + "execution_count": null, |
| 399 | + "metadata": { |
| 400 | + "tags": [] |
| 401 | + }, |
| 402 | + "outputs": [], |
| 403 | + "source": [ |
| 404 | + "rgb_during('2023-01-01')" |
| 405 | + ] |
| 406 | + }, |
369 | 407 | { |
370 | 408 | "cell_type": "markdown", |
371 | 409 | "metadata": {}, |
|
380 | 418 | "outputs": [], |
381 | 419 | "source": [ |
382 | 420 | "winter = '2022-01-15'\n", |
383 | | - "spring = '2022-04-15'\n", |
| 421 | + "spring = '2022-04-30'\n", |
384 | 422 | "summer = '2022-08-01'\n", |
385 | 423 | "fall = '2022-09-15'\n", |
386 | 424 | "\n", |
|
0 commit comments