07_Mapeo de severidad con Python y Sentinel-2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src='./IMG/head_b.png' alt='Logo Head' align='center' width='100%'></img>\n",
    "<br>\n",
    "<img src='./IMG/header_1.png' alt='Logo UNSPIDER' align='left' width='50%'></img>\n",
    "<br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a href=\"./05_MAPEO_INTERACTIVO_USANDO_API_PYTHON_DE_GEE_GEEMAP.ipynb\">**<< 05 - MAPA INTERACTIVO USANDO API Python de geemap**</a><span style=\"float:right;\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 7. Mapeo de severidad con Python y Sentinel-2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "http://www.un-spider.org/advisory-support/recommended-practices/recommended-practice-burn-severity/Step-by-Step/python-sentinel2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "El objetivo de esta práctica recomendada es generar un mapa de severidad de quemaduras para la evaluación de las áreas afectadas por los incendios forestales utilizando la Relación de quema normalizada diferenciada (dNBR), que es la diferencia entre la NBR posterior al incendio y la NBR previa al incendio."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Importar los paquetes requeridos"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "ename": "ModuleNotFoundError",
     "evalue": "No module named 'osr'",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mModuleNotFoundError\u001b[0m                       Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-9-112b391ef9b3>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[1;32mimport\u001b[0m \u001b[0mosr\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      2\u001b[0m \u001b[1;32mimport\u001b[0m \u001b[0mogr\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      3\u001b[0m \u001b[1;32mimport\u001b[0m \u001b[0mgdal\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      4\u001b[0m \u001b[1;32mimport\u001b[0m \u001b[0mnumpy\u001b[0m \u001b[1;32mas\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      5\u001b[0m \u001b[1;32mimport\u001b[0m \u001b[0mmatplotlib\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mModuleNotFoundError\u001b[0m: No module named 'osr'"
     ]
    }
   ],
   "source": [
    "import osr\n",
    "import ogr\n",
    "import gdal\n",
    "import numpy as np\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import rasterio\n",
    "import glob\n",
    "from rasterio.plot import show\n",
    "from rasterio.mask import mask\n",
    "from shapely.geometry import mapping\n",
    "import geopandas as gpd\n",
    "import math"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Definiciones de funciones"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def read_band_image(band, path):\n",
    "    \"\"\"\n",
    "    This function takes as input the Sentinel-2 band name and the path of the \n",
    "    folder that the images are stored, reads the image and returns the data as\n",
    "    an array\n",
    "    input:   band           string            Sentinel-2 band name\n",
    "             path           string            path of the folder\n",
    "    output:  data           array (n x m)     array of the band image\n",
    "             spatialRef     string            projection \n",
    "             geoTransform   tuple             affine transformation coefficients\n",
    "             targetprj                        spatial reference\n",
    "    \"\"\"\n",
    "    a = path+'*B'+band+'*.jp2'\n",
    "    img = gdal.Open(glob.glob(a)[0])\n",
    "    data = np.array(img.GetRasterBand(1).ReadAsArray())\n",
    "    spatialRef = img.GetProjection()\n",
    "    geoTransform = img.GetGeoTransform()\n",
    "    targetprj = osr.SpatialReference(wkt = img.GetProjection())\n",
    "    return data, spatialRef, geoTransform, targetprj\n",
    "\n",
    "def nbr(band1, band2):\n",
    "    \"\"\"\n",
    "    This function takes an input the arrays of the bands from the read_band_image\n",
    "    function and returns the Normalized Burn ratio (NBR)\n",
    "    input:  band1   array (n x m)      array of first band image e.g B8A\n",
    "            band2   array (n x m)      array of second band image e.g. B12\n",
    "    output: nbr     array (n x m)      normalized burn ratio\n",
    "    \"\"\"\n",
    "    nbr = (band1 - band2) / (band1 + band2)\n",
    "    return nbr\n",
    "\n",
    "def dnbr(nbr1,nbr2):\n",
    "    \"\"\"\n",
    "    This function takes as input the pre- and post-fire NBR and returns the dNBR\n",
    "    input:  nbr1     array (n x m)       pre-fire NBR\n",
    "            nbr2     array (n x m)       post-fire NBR\n",
    "    output: dnbr     array (n x m)       dNBR\n",
    "    \"\"\"\n",
    "    dnbr = nbr1 - nbr2\n",
    "    return dnbr\n",
    "\n",
    "def reproject_shp_gdal(infile, outfile, targetprj):\n",
    "    \"\"\"\n",
    "    This function takes as input the input and output file names and the projection\n",
    "    in which the input file will be reprojected and reprojects the input file using\n",
    "    gdal\n",
    "    input:  infile     string      input filename\n",
    "            outfile    string      output filename\n",
    "            targetprj              projection (output of function read_band_image)\n",
    "    \"\"\"\n",
    "    ## reprojection with gdal \n",
    "    \n",
    "    driver = ogr.GetDriverByName(\"ESRI Shapefile\") \n",
    "    dataSource = driver.Open(infile, 1) # 0 means read-only. 1 means writeable.\n",
    "    layer = dataSource.GetLayer()\n",
    "    sourceprj = layer.GetSpatialRef()\n",
    "    transform = osr.CoordinateTransformation(sourceprj, targetprj)\n",
    "    \n",
    "    # Create the output shapefile\n",
    "    outDriver = ogr.GetDriverByName(\"Esri Shapefile\")\n",
    "    outDataSource = outDriver.CreateDataSource(outfile)\n",
    "    outlayer = outDataSource.CreateLayer('', targetprj, ogr.wkbPolygon)\n",
    "    outlayer.CreateField(ogr.FieldDefn('id', ogr.OFTInteger))\n",
    "    \n",
    "    #Iterate over Features\n",
    "    i = 0\n",
    "    for feature in layer:\n",
    "        transformed = feature.GetGeometryRef()\n",
    "        transformed.Transform(transform) #reproject geometry\n",
    "\n",
    "        geom = ogr.CreateGeometryFromWkb(transformed.ExportToWkb()) # create geometry from wkb (write geometry of reprojected geometry)\n",
    "        defn = outlayer.GetLayerDefn() #layer definition\n",
    "        feat = ogr.Feature(defn)  #create new feature\n",
    "        feat.SetField('id', i) #set id\n",
    "        feat.SetGeometry(geom) #set geometry\n",
    "        outlayer.CreateFeature(feat) \n",
    "        i += 1\n",
    "        feat = None\n",
    "        \n",
    "def array2raster(array, geoTransform, projection, filename):\n",
    "    \"\"\" \n",
    "    This function tarnsforms a numpy array to a geotiff projected raster\n",
    "    input:  array                       array (n x m)   input array\n",
    "            geoTransform                tuple           affine transformation coefficients\n",
    "            projection                  string          projection\n",
    "            filename                    string          output filename\n",
    "    output: dataset                                     gdal raster dataset\n",
    "            dataset.GetRasterBand(1)                    band object of dataset\n",
    "    \n",
    "    \"\"\"\n",
    "    pixels_x = array.shape[1]\n",
    "    pixels_y = array.shape[0]\n",
    "    \n",
    "    driver = gdal.GetDriverByName('GTiff')\n",
    "    dataset = driver.Create(\n",
    "        filename,\n",
    "        pixels_x,\n",
    "        pixels_y,\n",
    "        1,\n",
    "        gdal.GDT_Float64, )\n",
    "    dataset.SetGeoTransform(geoTransform)\n",
    "    dataset.SetProjection(projection)\n",
    "    dataset.GetRasterBand(1).WriteArray(array)\n",
    "    dataset.FlushCache()  # Write to disk.\n",
    "    return dataset, dataset.GetRasterBand(1)  #If you need to return, remenber to return  also the dataset because the band don`t live without dataset.\n",
    " \n",
    "def clip_raster(filename, shp):\n",
    "    \"\"\"\n",
    "    This function clips a raster based on a shapefile\n",
    "    input:  filename          string                input raster filename\n",
    "            shp               dataframe             input shapefile open with geopandas\n",
    "    output: clipped           array (1 x n x m)     clipped array \n",
    "            clipped_meta      dict                  metadata\n",
    "            cr_ext            tuple                 extent of clipped data\n",
    "            gt                tuple                 affine transformation coefficients\n",
    "    \"\"\"\n",
    "    inraster = rasterio.open(filename)\n",
    "    \n",
    "    extent_geojson = mapping(shp['geometry'][0])\n",
    "    clipped, crop_affine = mask(inraster, \n",
    "                                shapes=[extent_geojson], \n",
    "                                nodata = np.nan,\n",
    "                                crop=True)\n",
    "    clipped_meta = inraster.meta.copy()\n",
    "    # Update the metadata to have the new shape (x and y and affine information)\n",
    "    clipped_meta.update({\"driver\": \"GTiff\",\n",
    "                 \"height\": clipped.shape[0],\n",
    "                 \"width\": clipped.shape[1],\n",
    "                 \"transform\": crop_affine})\n",
    "    cr_ext = rasterio.transform.array_bounds(clipped_meta['height'], \n",
    "                                            clipped_meta['width'], \n",
    "                                            clipped_meta['transform'])\n",
    "    \n",
    "    # transform to gdal\n",
    "    gt = crop_affine.to_gdal()\n",
    "    \n",
    "    return clipped, clipped_meta, cr_ext, gt\n",
    "    \n",
    "def reclassify(array):\n",
    "    \"\"\"\n",
    "    This function reclassifies an array\n",
    "    input:  array           array (n x m)    input array\n",
    "    output: reclass         array (n x m)    reclassified array\n",
    "    \"\"\"\n",
    "    reclass = np.zeros((array.shape[0],array.shape[1]))\n",
    "    for i in range(0,array.shape[0]):\n",
    "        for j in range(0,array.shape[1]):\n",
    "            if math.isnan(array[i,j]):\n",
    "                reclass[i,j] = np.nan\n",
    "            elif array[i,j] < 0.1:\n",
    "                reclass[i,j] = 1\n",
    "            elif array[i,j] < 0.27:\n",
    "                 reclass[i,j] = 2\n",
    "            elif array[i,j] < 0.44:\n",
    "                 reclass[i,j] = 3\n",
    "            elif array[i,j] < 0.66:\n",
    "                 reclass[i,j] = 4\n",
    "            else:\n",
    "                reclass[i,j] = 5\n",
    "                \n",
    "    return reclass\n",
    "                     "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "En la siguiente celda, debe modificar las rutas de las imágenes y el archivo de forma"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'gdal' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-11-cd458b14e4f2>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m     17\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     18\u001b[0m \u001b[1;31m# Read the pre-fire band images\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 19\u001b[1;33m \u001b[1;33m(\u001b[0m\u001b[0mpre_fire_b8a\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mcrs\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mgeoTransform\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mtargetprj\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mread_band_image\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mband1\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mpath_prefire\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     20\u001b[0m \u001b[1;33m(\u001b[0m\u001b[0mpre_fire_b12\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mcrs\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mgeoTransform\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mtargetprj\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mread_band_image\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mband2\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mpath_prefire\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     21\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m<ipython-input-10-96695809c327>\u001b[0m in \u001b[0;36mread_band_image\u001b[1;34m(band, path)\u001b[0m\n\u001b[0;32m     12\u001b[0m     \"\"\"\n\u001b[0;32m     13\u001b[0m     \u001b[0ma\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mpath\u001b[0m\u001b[1;33m+\u001b[0m\u001b[1;34m'*B'\u001b[0m\u001b[1;33m+\u001b[0m\u001b[0mband\u001b[0m\u001b[1;33m+\u001b[0m\u001b[1;34m'*.jp2'\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 14\u001b[1;33m     \u001b[0mimg\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mgdal\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mOpen\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mglob\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mglob\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0ma\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     15\u001b[0m     \u001b[0mdata\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0marray\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mimg\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mGetRasterBand\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mReadAsArray\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     16\u001b[0m     \u001b[0mspatialRef\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mimg\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mGetProjection\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mNameError\u001b[0m: name 'gdal' is not defined"
     ]
    }
   ],
   "source": [
    "path_prefire = \"F:/Burn_Severity/S2A_MSIL2A_20161220T143742_N0204_R096_T18HYF_20161220T145131.SAFE/GRANULE/L2A_T18HYF_A007815_20161220T145131/IMG_DATA/R20m/\"\n",
    "path_postfire = \"F:/Burn_Severity/S2A_MSIL2A_20170218T143751_N0204_R096_T18HYF_20170218T145150.SAFE/GRANULE/L2A_T18HYF_A008673_20170218T145150/IMG_DATA/R20m/\"\n",
    "# Define shapefile\n",
    "infile_shp = \"F:/Burn_Severity/Empedrado_adm_boundary/Empedrado.shp\"\n",
    "# Define reprojected shapefile\n",
    "outfile_shp = \"F:/Burn_Severity/Empedrado_adm_boundary/projected.shp\"\n",
    "# name of the output dNBR raster file\n",
    "filename = \"F:/Burn_Severity/dNBR.tiff\"\n",
    "# name of clipped dNBR raster\n",
    "filename2 = \"F:/Burn_Severity/dNBR_clipped.tiff\"\n",
    "# path to save figure\n",
    "fname = \"F:/Burn_Severity/map.png\"\n",
    "\n",
    "# Sentinel-2 Bands used for NBR calculation \n",
    "band1 = '8A'\n",
    "band2 = '12'\n",
    "    \n",
    "# Read the pre-fire band images \n",
    "(pre_fire_b8a, crs, geoTransform, targetprj) = read_band_image(band1, path_prefire)\n",
    "(pre_fire_b12, crs, geoTransform, targetprj) = read_band_image(band2, path_prefire)\n",
    "    \n",
    "# Calculation of pre-fire NBR\n",
    "pre_fire_nbr = nbr(pre_fire_b8a.astype(int),pre_fire_b12.astype(int))\n",
    "\n",
    "# Read the post-fire band images\n",
    "(post_fire_b8a, crs, geoTransform, targetprj) = read_band_image(band1, path_postfire)\n",
    "(post_fire_b12, crs, geoTransform, targetprj) = read_band_image(band2, path_postfire)\n",
    "    \n",
    "# Calculation of post-fire NBR\n",
    "post_fire_nbr = nbr(post_fire_b8a.astype(int),post_fire_b12.astype(int))\n",
    "    \n",
    "# Calculation of dNBR\n",
    "DNBR = dnbr(pre_fire_nbr,post_fire_nbr)\n",
    "    \n",
    "# Reprojection of shapefile with gdal to match projection of Sentinel-2 images\n",
    "reproject_shp_gdal(infile_shp, outfile_shp, targetprj)\n",
    "    \n",
    "# Read the reprojected shapefile\n",
    "fire_boundary = gpd.read_file(outfile_shp)\n",
    "    \n",
    "# project dNBR to images projection\n",
    "dnbr_tif, dnbr_tifBand = array2raster(DNBR, geoTransform, crs, filename)\n",
    "    \n",
    "# clip raster dNBR file to Empedrado shapefile\n",
    "(clipped_dnbr, clipped_dnbr_meta, cr_extent, gt) = clip_raster(filename, fire_boundary)\n",
    "clipped_ds , clipped_ds_rasterband = array2raster(clipped_dnbr[0], gt, crs, filename2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Graficar los resultados"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'matplotlib' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-12-e89ca0dc5fba>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m      2\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      3\u001b[0m \u001b[1;31m#set colors for plotting and classes\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 4\u001b[1;33m \u001b[0mcmap\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mmatplotlib\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mcolors\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mListedColormap\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m'green'\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;34m'yellow'\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;34m'orange'\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;34m'red'\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;34m'purple'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      5\u001b[0m \u001b[0mcmap\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mset_over\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'purple'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      6\u001b[0m \u001b[0mcmap\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mset_under\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'white'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mNameError\u001b[0m: name 'matplotlib' is not defined"
     ]
    }
   ],
   "source": [
    "# plot \n",
    "\n",
    "#set colors for plotting and classes\n",
    "cmap = matplotlib.colors.ListedColormap(['green','yellow','orange','red','purple'])\n",
    "cmap.set_over('purple')\n",
    "cmap.set_under('white')\n",
    "bounds = [-0.5, 0.1, 0.27, 0.440, 0.660, 1.3]        \n",
    "norm = matplotlib.colors.BoundaryNorm(bounds, cmap.N)  \n",
    "\n",
    "fig, ax = plt.subplots(figsize=(10, 10), subplot_kw={'xticks': [], 'yticks': []})\n",
    "cax = ax.imshow(clipped_ds_rasterband.ReadAsArray(), cmap=cmap, norm = norm)\n",
    "plt.title('Burn Severity Map')\n",
    "cbar = fig.colorbar(cax, ax=ax, fraction=0.035, pad=0.04, ticks=[-0.2, 0.18, 0.35, 0.53, 1])\n",
    "cbar.ax.set_yticklabels(['Unburned', 'Low Severity', 'Moderate-low Severity', 'Moderate-high Severity', 'High Severity'])\n",
    "plt.show()\n",
    "plt.savefig(fname, bbox_inches=\"tight\") \n",
    "        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calcular área quemada"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'clipped_ds_rasterband' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-13-73cb5da14160>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m      1\u001b[0m \u001b[1;31m# calculate burnt area (pixel size 20m*20m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 2\u001b[1;33m \u001b[0mreclass\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mreclassify\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mclipped_ds_rasterband\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mReadAsArray\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      3\u001b[0m \u001b[0mk\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m[\u001b[0m\u001b[1;34m'Unburned hectares'\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;34m'Low severity hectares'\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;34m'Moderate-low severity hectares'\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;34m'Moderate-high severity hectares'\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;34m'High severity'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      4\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m6\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      5\u001b[0m     \u001b[0mx\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mreclass\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mreclass\u001b[0m \u001b[1;33m==\u001b[0m \u001b[0mi\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mNameError\u001b[0m: name 'clipped_ds_rasterband' is not defined"
     ]
    }
   ],
   "source": [
    "# calculate burnt area (pixel size 20m*20m)\n",
    "reclass = reclassify(clipped_ds_rasterband.ReadAsArray())\n",
    "k = ['Unburned hectares', 'Low severity hectares', 'Moderate-low severity hectares', 'Moderate-high severity hectares', 'High severity']\n",
    "for i in range(1,6):\n",
    "    x = reclass[reclass == i]\n",
    "    l= x.size*0.04\n",
    "    print(\"%s: %.2f\" % (k[i-1], l))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}