updated multidimensional notebook

03_multidim_analysis.ipynb

@@ -98,12 +98,25 @@
     "nirband = 'LC08_L1TP_227065_20200608_20200626_01_T1_B{}.TIF'.format(5)\n",
     "mtlfile = 'LC08_L1TP_227065_20200608_20200626_01_T1_{}.json'.format('MTL')\n",
     "\n",
+    "# Overviews are reduced resolution versions of your dataset that can speed up rendering when you don’t need \n",
+    "# full resolution. By precomputing the upsampled pixels, rendering can be significantly faster when zoomed out.\n",
+    "# More info here: https://rasterio.readthedocs.io/en/latest/topics/overviews.html\n",
+    "\n",
     "with rasterio.open(url+redband) as src:\n",
     "    profile = src.profile\n",
     "    oviews = src.overviews(1) # list of overviews from biggest to smallest\n",
     "    oview = oviews[1]  # Use second-highest resolution overview\n",
     "    print('Decimation factor= {}'.format(oview))\n",
-    "    red = src.read(1, out_shape=(1, int(src.height // oview), int(src.width // oview)))\n",
+    "    red = src.read(1, out_shape=(1, int(src.height // oview), int(src.width // oview)))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Let's look at what the red band image looks like\n",
     "\n",
     "plt.imshow(red)\n",
     "plt.colorbar()\n",
@@ -133,6 +146,10 @@
    "metadata": {},
    "outputs": [],
    "source": [
+    "nir_filename = url + nirband\n",
+    "red_filename = url + redband\n",
+    "mtl_filename = url + mtlfile\n",
+    "\n",
     "download_file(nir_filename, 'nir.tif')\n",
     "download_file(red_filename, 'red.tif')\n",
     "download_file(mtl_filename, 'meta.json')"
@@ -146,7 +163,7 @@
    "source": [
     "# Get the shape size for the red band image\n",
     "\n",
-    "red = rasterio.open(url+redband)\n",
+    "red = rasterio.open(red_filename)\n",
     "print(red.is_tiled)\n",
     "red.block_shapes"
    ]
@@ -157,8 +174,10 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "red = xr.open_rasterio(url+redband, chunks={'band': 1, 'x': 1024, 'y': 1024})\n",
-    "nir = xr.open_rasterio(url+nirband, chunks={'band': 1, 'x': 1024, 'y': 1024})\n",
+    "# Technically you do not need to download the files to read the data\n",
+    "\n",
+    "red = xr.open_rasterio(red_filename, chunks={'band': 1, 'x': 1024, 'y': 1024})\n",
+    "nir = xr.open_rasterio(nir_filename, chunks={'band': 1, 'x': 1024, 'y': 1024})\n",
     "red"
    ]
   },
@@ -191,7 +210,9 @@
    "source": [
     "The Landsat Level 1 images are delivered in a quantized format. This has to be converted to top-of-atmosphere reflectance using the provided metadata.\n",
     "\n",
-    "First we define convenience functions to load the rescaling factors and transform a dataset. The red band is band 4 and near infrared is band 5."
+    "First we define convenience functions to load the rescaling factors and transform a dataset. The red band is band 4 and near infrared is band 5.\n",
+    "\n",
+    "We will also introduce here the concept of Dask, a flexible library for parallel computing in Python"
    ]
   },
   {
@@ -201,6 +222,18 @@
     "N𝐷𝑉𝐼=𝑁𝐼𝑅−𝑅𝑒𝑑 / 𝑁𝐼𝑅+𝑅𝑒𝑑"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import dask\n",
+    "from dask.distributed import Client\n",
+    "client = Client(processes=False)\n",
+    "client"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -237,42 +270,64 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "ndvi = (nir - red) / (nir + red)\n",
-    "ndvi2d = ndvi.squeeze()"
+    "red_toa = calculate_reflectance(red, band_number=4)\n",
+    "nir_toa = calculate_reflectance(nir, band_number=5)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "print(red_toa.variable.data)"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
    "metadata": {},
    "outputs": [],
-   "source": []
+   "source": [
+    "red_max, red_min, red_mean = dask.compute(\n",
+    "    red_toa.max(dim=['x', 'y']),\n",
+    "    red_toa.min(dim=['x', 'y']),\n",
+    "    red_toa.mean(dim=['x', 'y'])\n",
+    ")\n",
+    "print(red_max.item())\n",
+    "print(red_min.item())\n",
+    "print(red_mean.item())"
+   ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
    "metadata": {},
    "outputs": [],
-   "source": []
+   "source": [
+    "ndvi = (nir_toa - red_toa) / (nir_toa + red_toa)"
+   ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
-    "plt.figure()\n",
-    "im = ndvi2d.compute().plot.imshow(cmap='BrBG', vmin=-0.5, vmax=1)\n",
-    "plt.axis('equal')\n",
-    "plt.show()"
+    "ndvi2d = ndvi.squeeze()"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
    "metadata": {},
    "outputs": [],
-   "source": []
+   "source": [
+    "plt.figure()\n",
+    "im = ndvi2d.compute().plot.imshow(cmap='BrBG', vmin=-0.5, vmax=1)\n",
+    "plt.axis('equal')\n",
+    "plt.show()"
+   ]
   }
  ],
  "metadata": {