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

examples/bk/download_images

+#!/bin/python
+# Copyright (C) 2018 Adan Salazar <asalazargaribay@gmail.com>
+#
+#
+# This file is part of GeoSentinel
+#
+#
+# GeoSentinel is free software  you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# GeoSentinel is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with StereoVision.  If not, see <http://www.gnu.org/licenses/>.
+
+"""
+Example of dowloading Sentinel images from the Copernicus Open Access Hub < https://scihub.copernicus.eu/dhus/#/home >
+"""
+from sentinelsat.sentinel import read_geojson, geojson_to_wkt
+from argparse import ArgumentParser
+from geosentinel.ui_utils import (download_sentinel_mages)
+from geosentinel.arguments import SENTINEL_ARGUMENTS
+
+def main():
+    """   
+    Download all Sentinel images from the Copernicus Hub and save them in the folder specified by the user.
+    
+    Please run: python dowload_images  --help   before usage
+    
+    First, parse arguments provided by the user. Then choose the images that are in the
+    dates and polygon provided by the user, then remove the images that do not overlap more than 30 %
+    of the input polygon and the polygon of the image. Finally, download the remaining images.
+    
+    """
+    parser = ArgumentParser(description="Download Sentinel images from the Copernicus Hub",
+                             parents=[SENTINEL_ARGUMENTS])
+       
+     
+    args = parser.parse_args()
+    
+    path_json=args.input_file
+    
+    args.footprint = geojson_to_wkt(read_geojson(path_json))
+    
+    download_sentinel_mages(args)
+
+
+if __name__ == "__main__":
+    main()

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")
+

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/bk/JP2ToJPGSubset.py

+# -*- coding: utf-8 -*-
+# @package JP2ToJPGSubSet Sentinel Files
+"""
+Example:
+
+$ python -m doctest -v JP2ToJPGSubSet.py
+"""
+
+import os,sys,numpy as np
+from tools import (tools)
+from osgeo import gdal,ogr,osr
+
+
+class JP2ToJPG_SubSet(object):
+    """ Class for Sentinel satellites configuration
+    Test Case
+    >>> Files = tools("all") #IW_SLC,IW_GRDH,MSIL1C
+    >>> dir="make/"
+    >>> platform='Sentinel-2'
+    >>> list_files=Files.findFilesSentinel(dir,platform,2)#4 for resample_subset.dim
+    >>> print 'There is/are %d ' % len(list_files) + platform + ' file(s)'
+    >>> lat1=20.795
+    >>> lon1=-103.714
+    >>> lat2=20.165
+    >>> lon2=-103.086
+    >>> TypeMerge="RGB"
+    >>> Convert=JP2ToJPG_SubSet(TypeMerge)
+    >>> for image in list_files:
+    ... file = dir + "/"+ image
+    ... products = Files.search_band_combination(file,TypeMerge)
+    ... print"Applying Convert JP2 To JPEG and SubSet... ", image
+    ... Convert.Convert_SubSet(lat1, lon1, lon2, lat2, products)
+    """
+
+    def __init__(self,TypeMerge):
+        """The constructor Initialize Sentinel Data.
+
+        Args:
+        self: The object pointer.
+        TypeMerge (str): RGB
+
+        Returns:
+        pointer: The object pointer.
+        """
+        self.ext=".jpg"
+        self.SubSet="_SubSet"
+        self.Type_Merge="_"+TypeMerge
+
+    def Convert_SubSet(self, lat1, lon1, lon2, lat2, products):
+        """Index Files Sentinel
+
+         Args:
+             self (pointer): The object pointer.
+             lat1 (float): World Geodetic System 1984.
+             lon1 (float): World Geodetic System 1984.
+             lon2 (float): World Geodetic System 1984.
+             lat2 (float): World Geodetic System 1984.
+             files (str): local path to image.
+
+         """
+        for i in products:
+            outDS=i[:-4]+self.SubSet+self.ext
+            translate = 'gdal_translate -projwin %s %s %s %s -projwin_srs WGS84 -ot Byte -scale 0 4096 0 255 -of JPEG %s %s' % (lon1, lat1, lon2, lat2, i, outDS)
+            os.system(translate)
\ No newline at end of file

geosentinel/ndvi.py

+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+#import cv2
+import sys
+import numpy as np
+import gdal
+from osgeo import osr
+
+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)
+	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()
+	sr = osr.SpatialReference()
+	sr.ImportFromWkt(wkt)
+
+	output.SetProjection(sr.ExportToWkt())
+	output = None
+
+
+def main(argv):
+	geoint_pr_ndvi(argv[1], argv[2], argv[3])
+
+if __name__ == "__main__":
+	main(sys.argv)

geosentinel/polygonToBox.py

+#!/usr/bin/python
+# -*- coding: utf-8 -*-
+import sys
+from osgeo import ogr
+
+
+def getWKTPolygonBoundingBox(polygon):
+	'''!
+		Get the bounding box of a WKT polygon.
+		
+		@param polygon: WKT polygon
+		@return bounding box string in format ulx uly lrx lry
+
+		>>> g = getWKTPolygonBoundingBox("POLYGON((-89.96360136135893 20.754157792347172,-89.24078398227726 20.754157792347172,-89.24078398227726 21.295577631081912,-89.96360136135893 21.295577631081912,-89.96360136135893 20.754157792347172))")
+		>>> print (g)
+		-89.96360136135893 21.295577631081912 -89.24078398227726 20.754157792347172
+	'''
+
+	box = ogr.CreateGeometryFromWkt(polygon).GetEnvelope() #(minX, maxX, minY, maxY)
+
+	ulx = box[0]
+	uly = box[3]
+	lrx = box[1] 
+	lry = box[2]
+#	ulx uly lrx lry
+	return str(ulx) +" "+ str(uly) +" "+ str(lrx) +" "+ str(lry)
+
+def main(argv):
+#	print (argv[1])
+	print (getWKTPolygonBoundingBox(argv[1]))
+
+if __name__ == "__main__":
+	main(sys.argv)
+#g = getWKTPolygonBoundingBox("POLYGON((-89.96360136135893 20.754157792347172,-89.24078398227726 20.754157792347172,-89.24078398227726 21.295577631081912,-89.96360136135893 21.295577631081912,-89.96360136135893 20.754157792347172))")
+#print (g)

tools/productListToNDVI.sh

+#!/bin/sh
+
+BASEDIR=$1
+OUTDIR=$2
+WINDOW=$3
+count=1
+ext=.jpg
+BOX=$(python3 "$HOME/git/GeoSentinel/geosentinel/polygonToBox.py" "$3")
+echo $BOX
+
+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"08.jp2" -t |head -1)
+	VIR=$(ls $filepattern"04.jp2" -t |head -1)
+
+	echo "NIR-"$NIR
+	echo "VIR-"$VIR
+
+	if [ ! -d "ndvi" ]; then
+		mkdir ndvi
+	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"
+	fi
+
+        if [ ! -d "video_ndvi" ]; then
+                mkdir video_ndvi
+        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"
+#	set GDAL_PAM_ENABLED=YES
+	rm *.xml
+	cd ..
+
+done
+cd $OUTDIR"/video_ndvi"
+ffmpeg  -i %*.jpg -c:v libx264 -vf fps=10 -pix_fmt yuvj422p ndvi.mp4
+
+

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 
         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))
+
+        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
+
+
 	cd ..
 
 done
+
+
+cd out/video
+ffmpeg  -i %*.jpg -c:v libx264 -vf fps=10 -pix_fmt yuvj422p out.mp4
+