ndvi tiff

examples/ndviTest.py.save

-#!/usr/bin/python
-# -*- coding: utf-8 -*-
-import sys
-sys.path.append('../')
-from geosentinel import ndvi
-from datetime import date 
-
-
-ndvi.geoint_pr_ndvi(nir_fn, vir_fn, "out.jp2")
-

geosentinel/ndvi.py

@@ -17,7 +17,7 @@ def geoint_pr_ndvi(nir_fn, vir_fn, outfile):
 
 	rows, cols, geotransform = nir.RasterYSize, nir.RasterXSize, nir.GetGeoTransform()
 
-	# Read the input bands as numpy arrays.
+#	# Read the input bands as numpy arrays.
 	np_nir = nir.GetRasterBand(1).ReadAsArray(0,0,nir.RasterXSize, nir.RasterYSize) 
         np_red = red.GetRasterBand(1).ReadAsArray(0,0,red.RasterXSize, red.RasterYSize) 
 
@@ -30,33 +30,21 @@ def geoint_pr_ndvi(nir_fn, vir_fn, outfile):
 	numerator = np.subtract(np_nir_as32, np_red_as32)
 	denominator = np.add(np_nir_as32, np_red_as32)
 	result = np.divide(numerator, denominator)
+	result = np.multiply((result + 1), (2**7 - 1))
 
-	driverJP2 = gdal.GetDriverByName('JP2OpenJPEG')
 	driverTiff = gdal.GetDriverByName('GTiff')
-
-	ndvi_int8 = np.multiply((result + 1), (2**7 - 1))
-#        output = driverTiff.Create(outfile, cols, rows, 1, gdal.GDT_Byte)
-
-	output = driverJP2.CreateCopy(outfile, nir, strict=0)
-
-
-	raster = np.zeros((rows, cols), dtype=np.uint16)
-
-	output_band = output.GetRasterBand(1)
-        output_band.SetNoDataValue(-99)
-        output_band.WriteArray(raster)
-	#output.GetRasterBand(1).WriteArray(raster)
+#	driverJP2 = gdal.GetDriverByName('JP2OpenJPEG')
+#	driverJP2ECW = gdal.GetDriverByName('JP2ECW')
+#	raster = np.zeros((rows, cols), dtype=np.uint8)
+	output = driverTiff.Create(outfile+".tiff", cols, rows, 1, gdal.GDT_Byte)
+#	output = driverTiff.CreateCopy(outfile+".tiff", nir, strict=0)
+#	output = driverJP2.CreateCopy(outfile, nir, strict=0)
+        output.GetRasterBand(1).SetNoDataValue(-99)
+        output.GetRasterBand(1).WriteArray(result)
 	output.GetRasterBand(1).FlushCache()
-	output.GetRasterBand(1).WriteArray(raster)
-        output.FlushCache()
+	output.SetGeoTransform(geotransform)
 	output = None
-#	output.SetGeoTransform(geotransform)
-	output_band = None
-	nir = None
-
 
-	#del output
-	#ndvi_int8)
 
 def main(argv):
 	geoint_pr_ndvi(argv[1], argv[2], argv[3])

geosentinel/out_ndvi.jp2.aux.xml

-<PAMDataset>
-  <PAMRasterBand band="1">
-    <NoDataValue>-9.90000000000000E+01</NoDataValue>
-    <Metadata domain="IMAGE_STRUCTURE">
-      <MDI key="COMPRESSION">JPEG2000</MDI>
-      <MDI key="NBITS">15</MDI>
-    </Metadata>
-  </PAMRasterBand>
-</PAMDataset>