Merge branch 'master' of gitlab.geoint.mx:adan.salazar/GeoSentinel

geosentinel/ndvi.py → geosentinel/normalizedIndex.py

@@ -7,45 +7,40 @@ import numpy as np
 import gdal
 from osgeo import osr
 
-def geoint_pr_ndvi(nir_fn, vir_fn, outfile):
+def geoint_pr_normalizedIndex(b1_fn, b2_fn, outfile):
 	'''!
-		Normalized Difference Vegetation Index (NDVI)
-		@param nir_fr String: Near Infra Red band (S2_B08) Image file name
-		@param vir String: Red Band (S2_B04) Image file name
-		@param outfile String: NDVI image file name
+		Normalized Difference Index (ND*I)
+		@param b1_fn String: band #1 (S2_B*) Image file name
+		@param b2_fn String: Band #2 (S2_B*) Image file name
+		@param outfile String: processed image file name
 	'''
-	nir = gdal.Open(nir_fn)
-	red = gdal.Open(vir_fn)
+	b1 = gdal.Open(b1_fn)
+	b2 = gdal.Open(b2_fn)
 
-	rows, cols, geotransform = nir.RasterYSize, nir.RasterXSize, nir.GetGeoTransform()
+	rows, cols, geotransform = b1.RasterYSize, b1.RasterXSize, b1.GetGeoTransform()
 
-#	# 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) 
+	# Read the input bands as numpy arrays.
+	np_b1 = b1.GetRasterBand(1).ReadAsArray(0,0,b1.RasterXSize, b1.RasterYSize)
+	np_b2 = b2.GetRasterBand(1).ReadAsArray(0,0,b2.RasterXSize, b2.RasterYSize) 
 
 	# Convert the np arrays to 32-bit floating point to make sure division will occur properly.
-	np_nir_as32 = np_nir.astype(np.float32)
-	np_red_as32 = np_red.astype(np.float32)
+	np_b1_as32 = np_b1.astype(np.float32)
+	np_b2_as32 = np_b2.astype(np.float32)
 
 	# Calculate the NDVI formula.
 	np.seterr(divide='ignore', invalid='ignore')
-	numerator = np.subtract(np_nir_as32, np_red_as32)
-	denominator = np.add(np_nir_as32, np_red_as32)
+	numerator = np.subtract(np_b1_as32, np_b2_as32)
+	denominator = np.add(np_b1_as32, np_b2_as32)
 	result = np.divide(numerator, denominator)
 	result = np.multiply((result + 1), (2**7 - 1))
 
 	driverTiff = gdal.GetDriverByName('GTiff')
-#	driverJP2 = gdal.GetDriverByName('JP2OpenJPEG')
-#	driverJP2ECW = gdal.GetDriverByName('JP2ECW')
-#	raster = np.zeros((rows, cols), dtype=np.uint8)
 	output = driverTiff.Create(outfile, 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.SetGeoTransform(geotransform)
-	wkt = nir.GetProjection()
+	wkt = b1.GetProjection()
 	sr = osr.SpatialReference()
 	sr.ImportFromWkt(wkt)
 

tools/productListToNDVI.sh → tools/productListToIndex.sh

 #!/bin/sh
-
-BASEDIR=$1
-OUTDIR=$2
-WINDOW=$3
+# Get Normalized index (B1-B2)/(B1+B2) for all products in PRODCUTSDIR
+PRODCUTSDIR=$1	#Products directory
+JP2DIR=$2	#JP2 Directory
+WINDOW=$3	#Crop Window
+OUTPUTDIR=$4	#Process output
+$B1=$5		
+$B2=$6
 count=1
 ext=.jpg
-BOX=$(python3 "$HOME/git/GeoSentinel/geosentinel/polygonToBox.py" "$3")
+BOX=$(python3 "$HOME/git/GeoSentinel/geosentinel/polygonToBox.py" "$WINDOW")
 echo $BOX
 
-cd $BASEDIR
+cd $PRODCUTSDIR
 
 for f in $(find . -type f -name '*.zip')
 do 
 
 	filepattern=$(echo $f | cut -d"_" -f6)"_"$(echo $f | cut -d"_" -f3)"_B"
-	filename=$OUTDIR$filepattern
+	filename=$JP2DIR$filepattern
 	echo $filename
-	if [ ! -e $filename"04.jp2" ]; then
+	if [ ! -e $filename$B1".jp2" ]; then
 
-                unzip -n -j $f *B04.jp2 -d $OUTDIR 
+		unzip" -n -j $f *"$B2".jp2 -d " $JP2DIR 
 	fi
 
-	if [ ! -e $filename"08.jp2" ]; then
+	if [ ! -e $filename$B2".jp2" ]; then
 
-		unzip -n -j $f *B08.jp2 -d $OUTDIR 
+		unzip" -n -j $f *"$B1".jp2 -d " $JP2DIR 
 	fi
 
 	
-	cd $OUTDIR
-	NIR=$(ls $filepattern"08.jp2" -t |head -1)
-	VIR=$(ls $filepattern"04.jp2" -t |head -1)
+	cd $JP2DIR
+	$B1FN=$(ls $filepattern$B1".jp2" -t |head -1)
+	$B2FN=$(ls $filepattern$B1".jp2" -t |head -1)
 
-	echo "NIR-"$NIR
-	echo "VIR-"$VIR
+	echo "B1-"$B1FN
+	echo "B2-"$B2FN
 
-	if [ ! -d "ndvi" ]; then
-		mkdir ndvi
+	if [ ! -d $OUTPUTDIR ]; then
+		mkdir $OUTPUTDIR
 	fi
 
 	fileout1=$(echo $filepattern | cut -d"_" -f1)
         fileout2=$(echo $filepattern | cut -d"_" -f2)
 
-	if [ ! -e ndvi/$fileout2"_"$fileout1".tiff" ]; then
-		python3 $HOME/git/GeoSentinel/geosentinel/ndvi.py $NIR $VIR ndvi/$fileout2"_"$fileout1".tiff"
+	if [ ! -e $OUTPUTDIR"/"$fileout2"_"$fileout1".tiff" ]; then
+		python3 $HOME/git/GeoSentinel/geosentinel/normalizedIndex.py $B1FN $B2FN $OUTPUTDIR"/"$fileout2"_"$fileout1".tiff"
 	fi
 
-        if [ ! -d "video_ndvi" ]; then
-                mkdir video_ndvi
+        if [ ! -d "video" ]; then
+                mkdir video
         fi
 
 #	set GDAL_PAM_ENABLED=NO
-	gdal_translate -projwin $BOX -projwin_srs WGS84 -ot Byte -scale 0 255 0 255 -of JPEG ndvi/$fileout2"_"$fileout1".tiff" video_ndvi/$fileout2"_"$fileout1".jpg"
+	gdal_translate -projwin $BOX -projwin_srs WGS84 -ot Byte -scale 0 255 0 255 -of JPEG $OUTPUTDIR"/"$fileout2"_"$fileout1".tiff" video/$fileout2"_"$fileout1".jpg"
 #	set GDAL_PAM_ENABLED=YES
 	rm *.xml
-	cd ..
+	
+	cd $PRODCUTSDIR
 
 done
-cd $OUTDIR"/video_ndvi"
+cd $JP2DIR"/video_ndvi"
 ffmpeg  -i %*.jpg -c:v libx264 -vf fps=10 -pix_fmt yuvj422p ndvi.mp4