alfin

examples/downloadTest.py

@@ -9,4 +9,11 @@ area = "POLYGON((-89.79030210118333 21.122657323983717,-89.77308220413153 21.122
 footprint = "Intersects(POLYGON((-89.79030210118333 21.122657323983717,-89.77308220413153 21.122657323983717,-89.77308220413153 21.140540053466538,-89.79030210118333 21.140540053466538,-89.79030210118333 21.122657323983717)))"
 products = sentinel.getProducts(area, ('20150101', '20180517'), {"platformname":"Sentinel-2", "cloudcoverpercentage":"[0 TO 10]"})
 print len(products)
+for p in products:
+	print products[p]['filename']
+print len(products)
+products=sentinel.filterProducts(products)
+for p in products:
+        print products[p]['filename']
+print len(products)
 #sentinel.downloadProducts(products,dir)

examples/ndviTest.py

+#!/usr/bin/python
+# -*- coding: utf-8 -*-
+import sys
+sys.path.append('../')
+from geosentinel import ndvi
+from datetime import date 
+
+
+ndvi.geoint_pr_ndvi("/home/geoint/sentinelImages/out/T16QBJ_20180515T161829_B02.jp2", "/home/geoint/sentinelImages/out/T16QBJ_20180515T161829_B04.jp2", "out.jp2")
+

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/APISentinel.py

@@ -79,11 +79,18 @@ class APISentinel(object):
             dir (str): destination directory.
         """
         os.chdir(dir)
-        self.api.download_all(products)
+#        self.api.download_all(products)
+	for p in products:
+		self.api.download(p)
 
 
-#    def filterProducts(self, products_list, user_footprint, threshold):
-#        products_down = OrderedDict()
+    def filterProducts(self, productList):
+#, user_footprint, threshold):
+
+        products = productList.copy()
+	for p in productList:
+		if productList[p]['filename'].find("OPER_PRD") != -1:
+			del products[p]
 #        products_down = products_list.copy()
 #        products_df = self.api.to_dataframe(products_list)
 
@@ -104,6 +111,6 @@ class APISentinel(object):
 #                # Deleting element
 #                del products_down[products_df.uuid[i]]
 
-#        return products_down
+        return products
 
 

geosentinel/ndvi.py

+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+#import cv2
+import sys
+import numpy as np
+import gdal
+def geoint_pr_ndvi(nir_fn, vir_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
+	'''
+	nir = gdal.Open(nir_fn)
+	red = gdal.Open(vir_fn)
+
+	rows, cols, geotransform = nir.RasterYSize, nir.RasterXSize, nir.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) 
+
+	# 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)
+
+	# 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)
+	result = np.divide(numerator, denominator)
+
+	driverJP2 = gdal.GetDriverByName('JP2OpenJPEG')
+	ndvi_int8 = np.multiply((result + 1), (2**7 - 1))
+	output = driverJP2.CreateCopy(outfile, nir, strict=0)
+
+        output_band = output.GetRasterBand(1)
+        output_band.SetNoDataValue(-99)
+        output_band.WriteArray(ndvi_int8)
+
+def main(argv):
+	geoint_pr_ndvi(argv[1], argv[2], argv[3])
+
+if __name__ == "__main__":
+	main(sys.argv)

tools/productListToNDVI.sh

+#!/bin/sh
+
+BASEDIR=$1
+OUTDIR=$2
+WINDOW=$3
+count=1
+ext=.jpg
+cd $BASEDIR
+for f in $(find . -type f -name '*.zip')
+do 
+
+	filepattern=$(echo $f | cut -d"_" -f6)"_"$(echo $f | cut -d"_" -f3)"_B"
+	filename=$OUTDIR$filepattern
+	echo $filename
+	if [ ! -e $filename"04.jp2" ]; then
+
+                unzip -n -j $f *B04.jp2 -d $OUTDIR 
+	fi
+
+	if [ ! -e $filename"08.jp2" ]; then
+
+		unzip -n -j $f *B08.jp2 -d $OUTDIR 
+	fi
+
+	
+	cd $OUTDIR
+	NIR=$(ls $filepattern"04.jp2" -t |head -1)
+	VIR=$(ls $filepattern"08.jp2" -t |head -1)
+
+	echo "NIR-"$NIR
+	echo "VIR-"$VIR
+
+
+
+#	if [  ! -e $filepattern"04.jpg" ]; then
+#		gdal_translate -projwin $3 -projwin_srs WGS84 -ot Byte -scale 0 4096 0 255 -of JPEG $NIR ndvi/$filepattern"04.jpg"
+#	fi
+
+#	if [  ! -e $filepattern"08jpg" ]; then
+#		gdal_translate -projwin $3 -projwin_srs WGS84 -ot Byte -scale 0 4096 0 255 -of JPEG $VIR ndvi/$filepattern"08.jpg"
+#	fi
+
+
+	if [ ! -d "ndvi" ]; then
+		mkdir ndvi
+	fi
+
+	fileout1=$(echo $filepattern | cut -d"_" -f1)
+        fileout2=$(echo $filepattern | cut -d"_" -f2)
+
+	python /home/geoint/git/GeoSentinel/geosentinel/ndvi.py $filepattern"04.jp2" $filepattern"08.jp2" ndvi/$fileout2"_"$fileout1".jpg"
+#	convert $filepattern"04.jpg" $filepattern"03.jpg" $filepattern"02.jpg" -resize 640x480\! -combine video/$fileout2"_"$fileout1".jpg"
+
+        if [ ! -d "video" ]; then
+                mkdir video
+        fi
+
+	gdal_translate -projwin $3 -projwin_srs WGS84 -ot Byte -scale 0 4096 0 255 -of JPEG ndvi/$fileout2"_"$fileout1".jpg" video/$fileout2"_"$fileout1".jpg"
+	
+	cd ..
+
+done

tools/sentinelToRGBJPG.sh

@@ -27,39 +27,50 @@ do
 	echo "red-"$red
 	echo "gre-"$green
 	echo "blu-"$blue
-#    m=$(pwd)
-#    red=$m/$red
-#    green=$m/$green
-#    blue=$m/$blue
 
-	if [  ! -e $filepattern"04.jpg" ]; then
-		jp2ToJPGandStretch $red
-	fi
+#	if [  ! -e $filepattern"04.jpg" ]; then
+#		gdal_translate -projwin $3 -projwin_srs WGS84 -ot Byte -scale 0 4096 0 255 -of JPEG $red $filepattern"04.jpg"
+#	fi
+
+#	if [  ! -e $filepattern"03.jpg" ]; then
+#		gdal_translate -projwin $3 -projwin_srs WGS84 -ot Byte -scale 0 4096 0 255 -of JPEG $green $filepattern"03.jpg"
+#	fi
 
-	if [  ! -e $filepattern"03.jpg" ]; then
-		jp2ToJPGandStretch $green
+#        if [  ! -e $filepattern"02.jpg" ]; then
+#		gdal_translate -projwin $3 -projwin_srs WGS84 -ot Byte -scale 0 4096 0 255 -of JPEG $blue $filepattern"02.jpg"
+#	fi
+
+	if [ ! -d "rgb" ]; then
+		mkdir rgb
 	fi
 
-        if [  ! -e $filepattern"02.jpg" ]; then
-		jp2ToJPGandStretch $blue
+	fileout1=$(echo $filepattern | cut -d"_" -f1)
+        fileout2=$(echo $filepattern | cut -d"_" -f2)
+
+#	convert $filepattern"04.jpg" $filepattern"03.jpg" $filepattern"02.jpg" -resize 640x480\! -combine video/$fileout2"_"$fileout1".jpg"
+
+	if [  ! -e rgb/$fileout2"_"$fileout1".jp2" ]; then
+		gdal_merge.py -separate -co PHOTOMETRIC=RGB -o rgb/$fileout2"_"$fileout1".jp2" $filepattern"04.jp2" $filepattern"03.jp2" $filepattern"02.jp2"
 	fi
-#    red=$(ls *B04.jpg -t |head -1)
-#    green=$(ls *B03.jpg -t |head -1)
-#    blue=$(ls *B02.jpg -t |head -1)
-#    m=$(pwd)
-#    red=$m/$red
-#    green=$m/$green
-#    blue=$m/$blue
-	if [ ! -d "video" ]; then
-		mkdir video
+
+
+
+        if [ ! -d "video" ]; then
+                mkdir video 
+        fi
+
+        if [  ! -e video/$fileout2"_"$fileout1".jpg" ]; then
+
+	        gdal_translate -projwin $3 -projwin_srs WGS84 -ot Byte -scale 0 4096 0 255 -of JPEG rgb/$fileout2"_"$fileout1".jp2" video/$fileout2"_"$fileout1".jpg"
+	        convert video/$fileout2"_"$fileout1".jpg" -resize 640x480\!  video/$fileout2"_"$fileout1".jpg"
 	fi
-	fileout=$(echo $filepattern | cut -d"_" -f2)
-	convert $filepattern"04.jpg" $filepattern"02.jpg" $filepattern"02.jpg"  -combine video/$fileout".jpg"
-#    rm *.jp2
-#    rm *B04.jpg
-#    rm *B03.jpg
-#    rm *B02.jpg
-#    count=$((count+1))
+
+
 	cd ..
 
 done
+
+
+cd out/video
+ffmpeg  -i %*.jpg -c:v libx264 -vf fps=10 -pix_fmt yuvj422p out.mp4
+