Skip to content

Instantly share code, notes, and snippets.

@merkato merkato/visVegInd.py
Last active Jun 26, 2017

Embed
What would you like to do?
Calculating ortophoto indices by rasterio and numpy
#!/usr/bin/python
# -*- coding: utf-8 -*-
import rasterio
import os, sys
import numpy as np
from funkcje import *
try:
file = sys.argv[1]
typ = sys.argv[2]
(fileRoot, fileExt) = os.path.splitext(file)
outFileName = fileRoot + "_" + typ + fileExt
except IndexError:
print 'Not all arguments passed'
usage()
sys.exit()
except NameError:
print 'Vegetation Index not supported'
usage()
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, 0.0)
return np.choose(mask, (-9999.0, np.subtract(green, np.multiply(0.39,red), np.multiply(0.61, blue))))
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)
#dst.write_colormap(1, )
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.