MSAVI.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import all libraries\n",
    "import rasterio\n",
    "#import rasterio.plot\n",
    "import numpy as np\n",
    "import numpy.ma as ma\n",
    "import matplotlib.pyplot as plt\n",
    "from glob import glob\n",
    "import os\n",
    "import sys"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Modified Soil Adjusted Vegetation Index (MSAVI) description\n",
    "### Reference: Qi, J., Chehbouni, A.L., Huete, A.R., Kerr, Y.H., Sorooshian, S., 1994. A modified soil adjusted vegetation index (MSAVI). Remote Sensing of Environment 48, 119–126.\n",
    "\n",
    "\n",
    "### General description\n",
    "-This code has been tested with Sentinel-2B MSI Level 2 satellite data. \n",
    "-Any Sentinel-2B and 2A MSI Level 2 data can be download from https://scihub.copernicus.eu/dhus/#/home  \n",
    "-Due to three different spectral resolution of Sentinel-2 data (10m, 20m, and 60m), a preprocessing as resmapling has been applied to make all bands with 10m resolution. The resampling procedure has been done by SNAP ESA software.\n",
    "-The resampled data has been subset by the SANP ESA software to reduce the data memory.\n",
    "-Resampling and subset can also be done by python API of SNAP called snappy. Details on use sanppy API can be found here: https://senbox.atlassian.net/wiki/spaces/SNAP/pages/19300362/How+to+use+the+SNAP+API+from+Python\n",
    "-Resampling with snappy is not straight forward, due to  \n",
    "-For understanding the band with specific radiometric resolution, see here https://sentinel.esa.int/web/sentinel/user-guides/sentinel-2-msi/resolutions/radiometric"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### open band 8 with central wavelength of approx. 832.9 nm or IR or Infrared band, \n",
    "### for Sentinel-2A Band 8 is 832.8 nm\n",
    "### for Sentinel-2B Band 8 is 832.9 nm\n",
    "ras1 = 'YOUR:/file_location/S2B_MSIL2A_20200708T102559_N0214_R108_*****_20200708T135942_resampled.data/B8.img'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### this part is not necessary, unless if someone would like to see this info\n",
    "with rasterio.open(ras1) as src1:\n",
    "    print(src1.profile) # print the profile of raster data\n",
    "    print(src1.count) # count the number of bands\n",
    "    print(src1.crs) # show the raster geocoordinate\n",
    "    print(src1.crs.wkt) #show the raster WKT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### read band as 2D numpy array\n",
    "with rasterio.open(ras1, 'r') as ras1_src:\n",
    "    print(ras1_src)\n",
    "    ras1_arr = ras1_src.read(1)\n",
    "    print(ras1_arr.shape)\n",
    "    print(ras1_arr)    \n",
    "### remove no data value from array\n",
    "ndv = 0\n",
    "ras1_arr_ndv = ma.masked_where(ras1_arr == ndv, ras1_arr)\n",
    "print(ras1_arr_ndv.shape)\n",
    "print(ras1_arr_ndv)\n",
    "print(np.min(ras1_arr_ndv))\n",
    "print(np.max(ras1_arr_ndv))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### open band 4 with central wavelength of approx. 560 nm or Red visible band,\n",
    "### for Sentinel-2A Band 4 is 664.6 nm\n",
    "### for Sentinel-2B Band 4 is 664.9 nm\n",
    "ras2 = 'YOUR:/file_location/S2B_MSIL2A_20200708T102559_N0214_R108_*****_20200708T135942_resampled.data/B4.img'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### this part is not necessary, unless if someone would like to see this info\n",
    "with rasterio.open(ras2) as src2:\n",
    "    print(src2.profile) # print the profile of raster data\n",
    "    print(src2.count) # number of bands\n",
    "    print(src2.crs) # show the raster geocoordinate\n",
    "    print(src2.crs.wkt) #show the raster WKT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### read band as 2D numpy array\n",
    "with rasterio.open(ras2, 'r') as ras2_src:\n",
    "    print(ras2_src)\n",
    "    ras2_arr = ras2_src.read(1)\n",
    "    print(ras2_arr.shape)\n",
    "    print(ras2_arr)\n",
    "### remove no data value from array\n",
    "ndv = 0\n",
    "ras2_arr_ndv = ma.masked_where(ras2_arr == ndv, ras2_arr)\n",
    "print(ras2_arr_ndv.shape)\n",
    "print(ras2_arr_ndv)\n",
    "print(np.min(ras2_arr_ndv))\n",
    "print(np.max(ras2_arr_ndv))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### MSAVI calculation\n",
    "print(\"Calculating MSAVI.....\")\n",
    "MSAVI = 0.5*((2*(ras1_arr_ndv*0.0001)+1)-(((2*(ras1_arr_ndv*0.0001)+1)-8*((ras1_arr_ndv*0.0001)-(ras2_arr_ndv*0.0001)))**0.5))\n",
    "# multiply 0.0001 to convert 16-bit int to 32-bit float\n",
    "# because the Sentinel-2 Bottom of Atmosphere reflectance values are quantified by 10000\n",
    "# information can be found in the metadata: Level-2A_DataStrip_ID>\n",
    "# Image_Data_Info>Radiometric_info>QUANTIFICATION_VALUES_LIST\n",
    "\n",
    "print(MSAVI) # to see the calculated 2D numpy array\n",
    "print(MSAVI.shape) # to see the row and column of calculated array\n",
    "print(np.min(MSAVI)) # minimum value\n",
    "print(np.max(MSAVI)) # maximum value\n",
    "### plot array\n",
    "plt.imshow(MSAVI)\n",
    "### plot colorbar\n",
    "plt.colorbar(shrink=0.5)\n",
    "plt.clim(np.min(MSAVI), np.max(MSAVI))\n",
    "### add tile\n",
    "plt.title('MSAVI {}'.format(MSAVI.shape))\n",
    "### add x and y label\n",
    "plt.xlabel('Column #')\n",
    "plt.ylabel('Row #')\n",
    "### make the output file directory\n",
    "fig_direc = 'YOUR:/output_file_directory/MSAVI'\n",
    "if not os.path.exists(fig_direc):\n",
    "    os.makedirs(fig_direc)\n",
    "### Save the figure as jpg format\n",
    "plt.savefig('YOUR:/output_file_directory/MSAVI/MSAVI.jpg', dpi=600) # dpi can be changed according to needs but 600 is high \n",
    "# enough and such good resolution\n",
    "\n",
    "### Save calculated index as GeoTiff file\n",
    "### Get metadata from raster 2\n",
    "with rasterio.open(ras2) as src:\n",
    "    meta = src.meta\n",
    "meta.update(dtype=rasterio.float32)\n",
    "\n",
    "### Create output folder and write output file in it as geotiff\n",
    "MSAVI_file = \"YOUR:/output_file_directory/MSAVI/MSAVI.tif\"\n",
    "os.makedirs(os.path.dirname(MSAVI_file), exist_ok=True)\n",
    "with rasterio.open(MSAVI_file, 'w', **meta) as dst:\n",
    "    dst.write(MSAVI.astype(rasterio.float32), 1)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}