Calculate SAVI using snappy

Import libraries



In [ ]:

    
import numpy
from snappy import Product, ProductData, ProductIO, ProductUtils, FlagCoding

Read product



In [ ]:

    
file = '/eodata/Sentinel-2/MSI/L1C/2017/09/14/S2A_MSIL1C_20170914T033531_N0205_R061_T51VUJ_20170914T033641.SAFE'

L=0.5

print("Reading...")
product = ProductIO.readProduct(file)
width = product.getSceneRasterWidth()
height = product.getSceneRasterHeight()
name = product.getName()
description = product.getDescription()
band_names = product.getBandNames()

print("Product:     %s, %s" % (name, description))
print("Raster size: %d x %d pixels" % (width, height))
print("Start time:  " + str(product.getStartTime()))
print("End time:    " + str(product.getEndTime()))
print("Bands:       %s" % (list(band_names)))


b7 = product.getBand('B4')
b10 = product.getBand('B8')

Calculate SAVI and write product



In [ ]:

    
saviProduct = Product('savi', 'savi', width, height)
saviBand = saviProduct.addBand('savi', ProductData.TYPE_FLOAT32)
saviFlagsBand = saviProduct.addBand('savi_flags', ProductData.TYPE_UINT8)
writer = ProductIO.getProductWriter('BEAM-DIMAP')

ProductUtils.copyGeoCoding(product, saviProduct)

saviFlagCoding = FlagCoding('savi_flags')
saviFlagCoding.addFlag("savi_LOW", 1, "savi below 0")
saviFlagCoding.addFlag("savi_HIGH", 2, "savi above 1")
group = saviProduct.getFlagCodingGroup()
group.add(saviFlagCoding)

saviFlagsBand.setSampleCoding(saviFlagCoding)

saviProduct.setProductWriter(writer)
saviProduct.writeHeader('snappy_savi_output.dim')

r7 = numpy.zeros(width, dtype=numpy.float32)
r10 = numpy.zeros(width, dtype=numpy.float32)

print("Writing...")

for y in range(height):
    print("processing line ", y, " of ", height)
    r7 = b7.readPixels(0, y, width, 1, r7)
    r10 = b10.readPixels(0, y, width, 1, r10)
    savi = (r10 - r7) / (r10 + r7+L)
    saviBand.writePixels(0, y, width, 1, savi)
    saviLow = savi < 0.0
    saviHigh = savi > 1.0
    saviFlags = numpy.array(saviLow + 2 * saviHigh, dtype=numpy.int32)
    saviFlagsBand.writePixels(0, y, width, 1, saviFlags)

saviProduct.closeIO()

print("Done.")