.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "example_gallery/raster_operations/ndvi.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_example_gallery_raster_operations_ndvi.py: .. _example ndvi: NDVI ==== Operation using two raster bands In this example we read two Sentinel 2 raster bands and combine them to calculate the Normalized Difference Vegetation Index (NDVI). Since these bands are already from the same Sentinel 2 acquisition, they are already 'aligned'. This means none of the bands have to be resampled before we combine them. .. tip :: If this example matches your use case, you may also be interested in the packages Rasterio or Rioxarray. This example is created to get the reader familiar with the synax of GridKit and to highlight some of the intricacies when working with this package. Also, a numpy warning is raised because we devide by zero for some cells. These turn into NaN-values. This is remedied in the script by replacing each NaN-value with 0 afterwards. .. GENERATED FROM PYTHON SOURCE LINES 24-33 .. code-block:: Python import matplotlib.pyplot as plt import numpy from gridkit import read_raster # Read in the bands required to determine the NDVI band_4 = read_raster("../../tests/data/20220708_S2_L2A_B04_raw.tiff").astype("int32") band_8 = read_raster("../../tests/data/20220708_S2_L2A_B08_raw.tiff").astype("int32") .. GENERATED FROM PYTHON SOURCE LINES 34-41 .. warning :: The source data is defined in `int16`, which is consequently also the numpy dtype after reading. The numbers that result from adding the two bands are too large to store in a dtype of `int16`. Numpy does not upcast automatically, so we need to be aware of the dtypes ourselves. In this case the dtype is increased to `int32`, by calling ``astype`` directly after reading. If this is not done, the script will run, but the resulting values will be incorrect. .. GENERATED FROM PYTHON SOURCE LINES 42-47 .. code-block:: Python # Determine the NDVI ndvi = (band_8 - band_4) / (band_8 + band_4) ndvi.data[~numpy.isfinite(ndvi)] = 0 # fill NaN values resulting from a division by 0 .. GENERATED FROM PYTHON SOURCE LINES 48-55 .. note :: If, for any cell, the denominator (band_8 + band_4) is zero, Numpy raises a warning. The return value for these cells is a ``numpy.nan``. They will show as white in the image. For a smoother looking visualization, these values are replaced with a zero in this example. This decision is of course use case dependent. .. GENERATED FROM PYTHON SOURCE LINES 56-67 .. code-block:: Python # Plot the result mpl_extent = (band_4.bounds[0], band_4.bounds[2], band_4.bounds[1], band_4.bounds[3]) fig, ax = plt.subplots(1) im_ndvi = ax.imshow(ndvi, cmap="RdYlGn", extent=mpl_extent, vmin=-1, vmax=1) fig.colorbar(im_ndvi, ax=ax, fraction=0.022, pad=0.01) ax.set_xlabel("lon") ax.set_ylabel("lat") ax.set_title(f"NDVI of scene in the alps \n EPSG:{ndvi.crs.to_epsg()}") plt.show() .. image-sg:: /example_gallery/raster_operations/images/sphx_glr_ndvi_001.png :alt: NDVI of scene in the alps EPSG:3857 :srcset: /example_gallery/raster_operations/images/sphx_glr_ndvi_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 3.920 seconds) .. _sphx_glr_download_example_gallery_raster_operations_ndvi.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: ndvi.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: ndvi.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: ndvi.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_