Skip to content

Instantly share code, notes, and snippets.

@merkato merkato/
Last active Jun 26, 2017

What would you like to do?
Calculating ortophoto indices by rasterio and numpy
# -*- coding: utf-8 -*-
import rasterio
import os, sys
import numpy as np
from funkcje import *
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'
except NameError:
print 'Vegetation Index not supported'
def calcNgrdi(red, green):
Normalized green red difference index
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(
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.
: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 =
profile = ds.profile
profile.update(dtype=rasterio.float32, count=1, nodata=-9999)
red = np.float32(
green = np.float32(
blue = np.float32(
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, '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.