Applied-GeoSolutions / gips

Geospatial Image Processing System
GNU General Public License v3.0
17 stars 5 forks source link

add tillage indices to MODIS driver #439

Closed ircwaves closed 6 years ago

ircwaves commented 6 years ago

For consistency with landsat and sentinel2, these should be added as separate products (not in a multi-band TIF):

Below is a snippet showing how to generate an NDVI GeoImage using gippy, but given that Rob is MODIS driver maintainer, I think it'd be fine to use the non-gippy style.

from time import time
t0 = time()
import gippy
from gippy.algorithms import Indices

gippy.Options.SetVerbose(5)
gippy.Options.SetNumCores(1)
gippy.Options.SetChunkSize(2)

GeoImage = gippy.GeoImage
vsipath_pat = (
    '/vsizip//var/gips-test-data/sentinel2/tiles/19TCH/2017030/'
    'S2A_MSIL1C_20170130T153551_N0204_R111_T19TCH_20170130T153712.zip/'
    'S2A_MSIL1C_20170130T153551_N0204_R111_T19TCH_20170130T153712.SAFE/'
    'GRANULE/L1C_T19TCH_A008402_20170130T153712/IMG_DATA/'
    'T19TCH_20170130T153551_B{}.jp2'
)

fns = [vsipath_pat.format(bstr) for bstr in ('04', '08')]
img = GeoImage(fns)

# GIPPY uses color-based band names to generically reference raster bands
for i, c in enumerate(['RED', 'NIR']):
    img.SetBandName(c, i + 1)

prodout = Indices(
    img,
    {'ndvi': '/var/gips-test-data/sentinel2/tiles/19TCH/2017030/IRC_test_ndvi.tif'}, dict(),)
t1 = time()
print('module timing: {:0.4f} seconds'.format(t1-t0))
ircwaves commented 6 years ago

Equations are encoded here: https://github.com/Applied-GeoSolutions/gippy/blob/4359158e15f6d3d35976cd78cde8beda2627ef16/GIP/GeoAlgorithms.cpp#L793-L802 :

...
                 } else if (prodname == "ndti") {
                    cimgout = (swir1-swir2).div(swir1+swir2);
                } else if (prodname == "crc") {
                    cimgout = (swir1-blue).div(swir2+blue);
                } else if (prodname == "crcm") {
                    cimgout = (swir1-green).div(swir2+green);
                } else if (prodname == "isti") {
                    cimgout = swir2.div(swir1);
                } else if (prodname == "sti") {
                    cimgout = swir1.div(swir2);
...