kopia lustrzana https://github.com/OpenDroneMap/ODM
Agricultural indexes and Sentera AGX710 renaming script
sbonaime
rodzic
086716d7d7
commit
f82e39463d
|
|
@ -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="<orthophoto.tif>",
|
||||
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="<outfile.tif>",
|
||||
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()
|
||||
|
||||
|
||||
|
|
@ -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="<orthophoto.tif>",
|
||||
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
|
||||
Ładowanie…
Reference in New Issue