Merge pull request #617 from merkato/patch-1

Initial commit for vegind.py

contrib/visveg/readme.md

@@ -0,0 +1,31 @@
+# Visible Vegetation Indexes
+
+This script produces a Vegetation Index raster from a RGB orthophoto (odm_orthophoto.tif in your project)
+
+## Requirements
+* rasterio (pip install rasterio)
+* numpy python package (included in ODM build)
+
+## Usage
+```
+vegind.py <orthophoto.tif> index
+
+positional arguments:
+  <orthophoto.tif>  The RGB orthophoto. Must be a GeoTiff.
+  index     Index identifier. Allowed values: ngrdi, tgi, vari
+```
+Output will be generated with index suffix in the same directory as input.
+
+## Examples
+
+`python vegind.py /path/to/odm_orthophoto.tif tgi`
+
+Orthophoto photo of Koniaków grass field and forest in QGIS: ![](http://imgur.com/K6x3nB2.jpg)
+The Triangular Greenness Index output in QGIS (with a spectral pseudocolor): ![](http://i.imgur.com/f9TzISU.jpg)
+Visible Atmospheric Resistant Index: ![](http://imgur.com/Y7BHzLs.jpg)
+Normalized green-red difference index: ![](http://imgur.com/v8cmaPS.jpg)
+
+## Bibliography
+
+1. Hunt, E. Raymond, et al. "A Visible Band Index for Remote Sensing Leaf Chlorophyll Content At the Canopy Scale." ITC journal 21(2013): 103-112. doi: 10.1016/j.jag.2012.07.020
+(https://doi.org/10.1016/j.jag.2012.07.020)

contrib/visveg/vegind.py

@@ -0,0 +1,94 @@
+#!/usr/bin/python
+# -*- coding: utf-8 -*-
+import rasterio, os, sys
+import numpy as np
+
+class bcolors:
+    OKBLUE = '\033[94m'
+    OKGREEN = '\033[92m'
+    WARNING = '\033[93m'
+    FAIL = '\033[91m'
+    ENDC = '\033[0m'
+    BOLD = '\033[1m'
+    UNDERLINE = '\033[4m'
+
+try:
+    file = sys.argv[1]
+    typ = sys.argv[2]
+    (fileRoot, fileExt) = os.path.splitext(file)
+    outFileName = fileRoot + "_" + typ + fileExt
+    isinstance(typ, ['vari', 'tgi', 'ngrdi'])
+except (TypeError, IndexError, NameError):
+    print bcolors.FAIL + 'Arguments messed up. Check arguments order and index name' + bcolors.ENDC
+    print 'Usage: ./vegind.py orto index'
+    print '       orto - filepath to RGB orthophoto'
+    print '       index - Vegetation Index'
+    print bcolors.OKGREEN + 'Available indexes: vari, ngrdi, tgi' + bcolors.ENDC
+    sys.exit()
+
+
+def calcNgrdi(red, green):
+    """
+    Normalized green red difference index
+    Tucker,C.J.,1979.
+    Red and photographic infrared linear combinations for monitoring vegetation.
+    Remote Sensing of Environment 8, 127–150
+    :param red: red visible channel
+    :param green: green visible channel
+    :return: ngrdi index array
+    """
+    mask = np.not_equal(np.add(red,green), 0.0)
+    return np.choose(mask, (-9999.0, np.true_divide(
+    np.subtract(green,red),
+    np.add(red,green))))
+
+def calcVari(red,green,blue):
+    """
+    Calculates Visible Atmospheric Resistant Index
+    Gitelson, A.A., Kaufman, Y.J., Stark, R., Rundquist, D., 2002.
+    Novel algorithms for remote estimation of vegetation fraction.
+    Remote Sensing of Environment 80, 76–87.
+    :param red: red visible channel
+    :param green: green visible channel
+    :param blue: blue visible channel
+    :return: vari index array, that will be saved to tiff
+    """
+    mask = np.not_equal(np.subtract(np.add(green,red),blue), 0.0)
+    return np.choose(mask, (-9999.0, np.true_divide(np.subtract(green,red),np.subtract(np.add(green,red),blue))))
+
+def calcTgi(red,green,blue):
+    """
+    Calculates Triangular Greenness Index
+    Hunt, E. Raymond Jr.; Doraiswamy, Paul C.; McMurtrey, James E.; Daughtry, Craig S.T.; Perry, Eileen M.; and Akhmedov, Bakhyt,
+    A visible band index for remote sensing leaf chlorophyll content at the canopy scale (2013).
+    Publications from USDA-ARS / UNL Faculty. Paper 1156.
+    http://digitalcommons.unl.edu/usdaarsfacpub/1156
+    :param red: red channel
+    :param green: green channel
+    :param blue: blue channel
+    :return: tgi index array, that will be saved to tiff
+    """
+    mask = np.not_equal(green-red+blue-255.0, 0.0)
+    return np.choose(mask, (-9999.0, np.subtract(green, np.multiply(0.39,red), np.multiply(0.61, blue))))
+
+try:
+    with rasterio.Env():
+        ds = rasterio.open(file)
+        profile = ds.profile
+        profile.update(dtype=rasterio.float32, count=1, nodata=-9999)
+        red = np.float32(ds.read(1))
+        green = np.float32(ds.read(2))
+        blue = np.float32(ds.read(3))
+        np.seterr(divide='ignore', invalid='ignore')
+        if typ == 'ngrdi':
+            indeks = calcNgrdi(red,green)
+        elif typ == 'vari':
+            indeks = calcVari(red, green, blue)
+        elif typ == 'tgi':
+            indeks = calcTgi(red, green, blue)
+
+        with rasterio.open(outFileName, 'w', **profile) as dst:
+            dst.write(indeks.astype(rasterio.float32), 1)
+except rasterio.errors.RasterioIOError:
+    print bcolors.FAIL + 'Orthophoto file not found or access denied' + bcolors.ENDC
+    sys.exit()