In [ ]:

    

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

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')
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()
#print(dir(group))
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.")


    

    




Reading...
Product:     S2A_MSIL1C_20170914T033531_N0205_R061_T51VUJ_20170914T033641, None
Raster size: 10980 x 10980 pixels
Start time:  14-SEP-2017 03:35:31.026000
End time:    14-SEP-2017 03:35:31.026000
Bands:       ['B1', 'B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'B8', 'B8A', 'B9', 'B10', 'B11', 'B12', 'view_zenith_mean', 'view_azimuth_mean', 'sun_zenith', 'sun_azimuth', 'view_zenith_B1', 'view_azimuth_B1', 'view_zenith_B2', 'view_azimuth_B2', 'view_zenith_B3', 'view_azimuth_B3', 'view_zenith_B4', 'view_azimuth_B4', 'view_zenith_B5', 'view_azimuth_B5', 'view_zenith_B6', 'view_azimuth_B6', 'view_zenith_B7', 'view_azimuth_B7', 'view_zenith_B8', 'view_azimuth_B8', 'view_zenith_B8A', 'view_azimuth_B8A', 'view_zenith_B9', 'view_azimuth_B9', 'view_zenith_B10', 'view_azimuth_B10', 'view_zenith_B11', 'view_azimuth_B11', 'view_zenith_B12', 'view_azimuth_B12']
Writing...
processing line  0  of  10980
processing line  1  of  10980
processing line  2  of  10980
processing line  3  of  10980
processing line  4  of  10980
processing line  5  of  10980
processing line  6  of  10980
processing line  7  of  10980
processing line  8  of  10980
processing line  9  of  10980
processing line  10  of  10980
processing line  11  of  10980
processing line  12  of  10980
processing line  13  of  10980
processing line  14  of  10980
processing line