diff --git a/contrib/ndvi/agricultural_indices.py b/contrib/ndvi/agricultural_indices.py new file mode 100755 index 00000000..0d3a40b6 --- /dev/null +++ b/contrib/ndvi/agricultural_indices.py @@ -0,0 +1,179 @@ +#!/usr/bin/env python3 +# A script to calculate agricultural indices +# NDVI - Normalized Difference Vegetation Index - (NIR−RED)/(NIR + RED) +# NDRE - Normalized Difference Red Edge - (NIR−RE)/(NIR + RE) +# GNDVI - Green NDVI - (NIR−GREEN)/(NIR + GREEN) +# GRVI - Green RVI - NIR/GREEN +# https://support.micasense.com/hc/en-us/articles/226531127-Creating-agricultural-indices-NDVI-NDRE-in-QGIS- +# requires python-gdal + + + +import numpy +import argparse +import os.path +try: + from osgeo import gdal + from osgeo import osr +except ImportError: + raise ImportError("You need to install python-gdal : \ + run `sudo apt-get install libgdal-dev` \ + # Check Gdal version with \ + gdal-config --version \ + #install correspondig gdal version with pip : \ + pip3 install GDAL==2.4.0") + + +def parse_args(): + argument_parser = argparse.ArgumentParser('Createa from a multispectral orthophoto \ +a Geotif with NDVI, NDRE, GNDVI and GRVI agricultural indices') + + argument_parser.add_argument("orthophoto", metavar="", + type=argparse.FileType('r'), + help="The CIR orthophoto. Must be a GeoTiff.") + argument_parser.add_argument("-red", type=int, + help="Red band number") + argument_parser.add_argument("-green", type=int, + help="Green band number") + argument_parser.add_argument("-blue", type=int, + help="Blue band number") + argument_parser.add_argument("-re", type=int, + help="RedEdge band number") + argument_parser.add_argument("-nir", type=int, + help="NIR band number") + argument_parser.add_argument("out", metavar="", + type=argparse.FileType('w'), + help="The output file.") + argument_parser.add_argument("--overwrite", "-o", + action='store_true', + default=False, + help="Will overwrite output file if it exists. ") + return argument_parser.parse_args() + + +def calc_ndvi(nir, red): + """ + Calculates the NDVI of an orthophoto using nir and red bands. + :param nir: An array containing the nir band + :param vis: An array containing the red band + :return: An array that will be exported as a tif + """ + + # Take the orthophoto and do nir - red / nir + red + # for each cell, calculate ndvi (masking out where divide by 0) + ndvi = numpy.empty(nir.shape, dtype=float) + mask = numpy.not_equal((nir + red), 0.0) + return numpy.choose(mask, (-1.0, numpy.true_divide(numpy.subtract(nir, red), numpy.add(nir, red)))) + +def calc_ndre(nir, re): + """ + Calculates the NDRE of an orthophoto using nir and re bands. + :param nir: An array containing the nir band + :param re: An array containing the rededge band + :return: An array that will be exported as a tif + """ + + # Take the orthophoto and do nir - re / nir + re + # for each cell, calculate ndre (masking out where divide by 0) + ndre = numpy.empty(nir.shape, dtype=float) + mask = numpy.not_equal((nir + re), 0.0) + return numpy.choose(mask, (-1.0, numpy.true_divide(numpy.subtract(nir, re), numpy.add(nir, re)))) + +def calc_gndvi(nir, green): + """ + Calculates the GNDVI of an orthophoto using nir and green bands. + :param nir: An array containing the nir band + :param green: An array containing the green band + :return: An array that will be exported as a tif + """ + + # Take the orthophoto and do nir - re / nir + re + # for each cell, calculate ndre (masking out where divide by 0) + gndvi = numpy.empty(nir.shape, dtype=float) + mask = numpy.not_equal((nir + green), 0.0) + return numpy.choose(mask, (-1.0, numpy.true_divide(numpy.subtract(nir, green), numpy.add(nir, green)))) + +def calc_grvi(nir, green): + """ + Calculates the GRVI of an orthophoto using nir and green bands. + :param nir: An array containing the nir band + :param green: An array containing the green band + :return: An array that will be exported as a tif + """ + + # Take the orthophoto and do nir - re / nir + re + # for each cell, calculate ndre (masking out where divide by 0) + grvi = numpy.empty(nir.shape, dtype=float) + mask = numpy.not_equal((nir + green), 0.0) + return numpy.choose(mask, (-1.0, numpy.true_divide(nir, green))) + + + +if __name__ == "__main__": + + # Supress/hide warning when dividing by zero + numpy.seterr(divide='ignore', invalid='ignore') + + rootdir = os.path.dirname(os.path.abspath(__file__)) + + # Parse args + args = parse_args() + + if not args.overwrite and os.path.isfile(os.path.join(rootdir, args.out.name)): + print("File exists, rename or use -o to overwrite.") + exit() + + # import raster + print("Reading file") + raster = gdal.Open(args.orthophoto.name) + orthophoto = raster.ReadAsArray() + + # parse out bands + print("Reading rasters") + red_matrix=orthophoto[args.red-1].astype(float) + green_matrix=orthophoto[args.green-1].astype(float) + blue_matrix=orthophoto[args.blue-1].astype(float) + re_matrix=orthophoto[args.re-1].astype(float) + nir_matrix=orthophoto[args.nir-1].astype(float) + + + outfile = args.out + + # NDVI + print("Computing NDVI") + #ndvi = calc_ndvi(nir_matrix, red_matrix) + ndvi = (nir_matrix.astype(float) - red_matrix.astype(float)) / (nir_matrix + red_matrix) + # NDRE + print("Computing NDRE") + #ndre = calc_ndre(nir_matrix, re_matrix) + ndre = (nir_matrix.astype(float) - re_matrix.astype(float)) / (nir_matrix + re_matrix) + + # GNDVI + print("Computing GNDVI") + #gndvi = calc_gndvi(nir_matrix, green_matrix) + gndvi = (nir_matrix.astype(float) - green_matrix.astype(float)) / (nir_matrix + green_matrix) + + # GRVI + print("Computing GRVI") + #grvi = calc_grvi(nir_matrix, green_matrix) + grvi = (nir_matrix.astype(float) / green_matrix) + + __import__("IPython").embed() + + print("Saving Files") + # export raster + + for name, matrix in zip(['ndvi', 'ndre', 'gndvi' ,'grvi'] ,[ndvi,ndre,gndvi,grvi] ): + print(name) + out_driver = gdal.GetDriverByName('GTiff')\ + .Create(name+'_'+outfile.name, int(ndvi.shape[1]), int(ndvi.shape[0]), 1, gdal.GDT_Float32) + outband = out_driver.GetRasterBand(1) + outband.SetDescription(name.capitalize()) + outband.WriteArray(matrix) + outcrs = osr.SpatialReference() + outcrs.ImportFromWkt(raster.GetProjectionRef()) + out_driver.SetProjection(outcrs.ExportToWkt()) + out_driver.SetGeoTransform(raster.GetGeoTransform()) + outband.FlushCache() + + diff --git a/contrib/ndvi/rename_sentera_agx710_multispectral_tif.py b/contrib/ndvi/rename_sentera_agx710_multispectral_tif.py new file mode 100644 index 00000000..d2077d16 --- /dev/null +++ b/contrib/ndvi/rename_sentera_agx710_multispectral_tif.py @@ -0,0 +1,73 @@ +#!/usr/bin/env python3 +# A script to rename. +# requires python-gdal + +import argparse +import sys +try: + from osgeo import gdal +except ImportError: + raise ImportError("You need to install python-gdal : \ + run `sudo apt-get install libgdal-dev` \ + # Check Gdal version with \ + gdal-config --version \ + #install correspondig gdal version with pip : \ + pip3 install GDAL==2.4.0") + +def parse_args(): + """ Parse arguments """ + argument_parser = argparse.ArgumentParser( + "A script that rename inplace Sentera AGX710 Geotiff orthophoto. ") + argument_parser.add_argument("orthophoto", metavar="", + type=argparse.FileType('r'), + help="The input orthophoto. Must be a GeoTiff.") + return argument_parser.parse_args() + + +def rename_sentera_agx710_layers(name): + """ Only rename Geotif built from Sentera AGX710 images with ODM """ + if raster.RasterCount != 7: + raise ImportError(F'File {name} does not have 7 layers as a regular\ + Geotif built from Sentera AGX710 images with ODM') + + if 'RedGreenBlue' in raster.GetRasterBand(1).GetDescription() and \ + 'RedEdgeGarbageNIR' in raster.GetRasterBand(2).GetDescription(): + + print("Sentera AGX710 Geotiff file has been detected.\ + Layers are name are :") + print("RedGreenBlue for Band 1\nRedEdgeGarbageNIR for Band 2\ + \nNone for Band 3\nNone for Band 4\nNone for Band 5\nNone for Band 6") + print("\nAfter renaming bands will be :") + print("Red for Band 1\nGreen for Band 2\nBlue for Band 3\n\ + RedEdge for Band 4\nGarbage for Band 5\nNIR for Band 6") + + answer = input( + "Are you sure you want to rename the layers of the input file ? [yes/no] ") + if answer =='yes': + raster.GetRasterBand(1).SetDescription('Red') + raster.GetRasterBand(2).SetDescription('Green') + raster.GetRasterBand(3).SetDescription('Blue') + raster.GetRasterBand(4).SetDescription('RedEdge') + raster.GetRasterBand(5).SetDescription('Garbage') + raster.GetRasterBand(6).SetDescription('NIR') + # raster.GetRasterBand(7).SetDescription('Alpha') + else: + print("No renaming") + else : + print(F'No need for band renaming in {name}') + sys.exit() + + +if __name__ == "__main__": + + # Parse args + args = parse_args() + + # import raster + raster = gdal.Open(args.orthophoto.name, gdal.GA_Update) + + # Rename layers + rename_sentera_agx710_layers(args.orthophoto.name) + + # de-reference the datasets, which triggers gdal to save + raster = None