update forom my laptop

geosentinel/APISentinel.py

@@ -22,7 +22,7 @@ import time
 class APISentinel(object):
 	""" Class for Sentinel satellites configuration
 	Test Case
-	>>> sentinel = APISentinel('asalazarg',  'geo135asg')
+	>>> sentinel = APISentinel('mario-chirinos',  'r4nc0r4u')
 	>>> products = sentinel.getProducts("POLYGON((-89.99450683593753 21.279137394108716,-89.13757324218751 21.2996106049456,-89.30236816406251 20.68418377935238,-90.0494384765625 20.715015145512098,-89.99450683593753 21.279137394108716))", ('20151219', date(2015,12,29)), {"platformname":"Sentinel-2"})
 	>>> print(len(products))
 	6
@@ -43,7 +43,8 @@ class APISentinel(object):
 		self.username = usernam
 		self.password = passw
 #		self.URL_Sentinel = 'https://scihub.copernicus.eu/dhus'
-		self.URL_Sentinel = 'https://scihub.copernicus.eu/apihub'
+#		self.URL_Sentinel = 'https://scihub.copernicus.eu/apihub'
+		self.URL_Sentinel = 'https://apihub.copernicus.eu/apihub/'
 	# self.plot_width = 7
 	# self.plot_height = 7
 	# self.intersection_th = 30

geosentinel/findProducts.py

@@ -16,7 +16,8 @@ from datetime import date
 from collections import OrderedDict
 
 productType = {"L1C":"S2MSI1C", "L2A":"S2MSI2A"}
-def findSentinelProducts(cfg, uinput = False):#main(argv):
+#===============================================================================
+def findSentinelProducts(cfg, uinput = False, limitto = 0):#main(argv):
 
 	sentinel = APISentinel.APISentinel(cfg['username'],  cfg['password'])
 	geometry = ogr.CreateGeometryFromWkt(cfg['wkt'])
@@ -37,7 +38,7 @@ def findSentinelProducts(cfg, uinput = False):#main(argv):
 
 	productList = OrderedDict()
 	for wkt in wktList:
-		tmpList	 =sentinel.getProducts(wkt, (cfg['startDate'], cfg['endDate']), {"producttype":productType[cfg["productLevel"]], "platformname":cfg['platform'], "cloudcoverpercentage":"[0 TO "+str(cfg['clouds'])+"]"}, "Contains")		
+		tmpList	 =sentinel.getProducts(wkt, (cfg['startDate'], cfg['endDate']), {"producttype":productType[cfg["productLevel"]], "platformname":cfg['platform'], "cloudcoverpercentage":"[0 TO "+str(cfg['clouds'])+"]"}, areaRelation="Contains")		
 
 		if len(productList)<=0:
 			print ("No products found with 'Contains' trying with 'Intersects'")
@@ -45,6 +46,21 @@ def findSentinelProducts(cfg, uinput = False):#main(argv):
 	
 		productList.update(tmpList)
 
+	if limitto>0:
+		tilesdic={}
+		for k, v in productList.items():
+#			print("productList v", v)
+			tile = v["title"][38:44]
+			date = v["title"][11:19]
+			if tile not in tilesdic:
+				tilesdic[tile]={}
+			tilesdic[tile][k]=v
+
+		for k,v in tilesdic.items():
+			tilesdic[k] = {i[0]:i[1] for i in sorted(v.items(), key=lambda x: x[1]["title"][11:19], reverse=True)}
+
+		productList = {list(tilesdic[k].keys())[0]: list(tilesdic[k].items())[0][1] for k in tilesdic.keys()}
+
 	fileNames = [productList[k]['filename'].replace("SAFE", "zip") for k in productList.keys() ]
 
 	downloadDir =  cfg["productsDir"]+cfg["productLevel"]+"/"

geosentinel/polygonToBox.py

@@ -35,11 +35,10 @@ def getWKTPolygonBoundingBox(polygon, wkt=False):
 		return str(ulx) +" "+ str(uly) +" "+ str(lrx) +" "+ str(lry)
 
 def main(argv):
-
-	if len(sys.argv) == 2:
+	if len(argv) == 2:
 		print (getWKTPolygonBoundingBox(argv[1]))
-	if len(sys.argv) == 3:
-		print (getWKTPolygonBoundingBox(argv[1], argv[2]))
+	if len(argv) == 3:
+		print (getWKTPolygonBoundingBox(argv[1], wkt=True))
 
 if __name__ == "__main__":
 	main(sys.argv)

geosentinel/rasterWkt.py

@@ -90,7 +90,7 @@ def rasterWkt(wkt, inputfile, outputfile):
 	polygonRaster = None
 	output = None
 	inputImage = None
-
+#-------------------------------------------------------------------------------
 def main(argv):
 	rasterWkt(argv[1], argv[2], argv[3])
 

tools/ExtractData.sh

 #!/bin/sh
+echo $PWD
 FILE=$1
 JP2DIR=$2	#JP2 Output Directory
 BYTILE=${3:-0}
 #S2A_MSIL2A_60170515T162341_N0605_R040_T16QBJ_60170515T163103.zip
 filepattern=$(echo $FILE | cut -d"_" -f6)"_"$(echo $FILE | cut -d"_" -f3)"_" #T16QBJ_60170515T162341_
 
-productType=$(echo $filepattern | cut -d"_" -f2)	#MSIL2A
+#productType=$(echo $filepattern | cut -d"_" -f2)	#MSIL2A
 tile=$(echo $filepattern | cut -d"_" -f1) 		#T16QBJ
 datetime=$(echo $filepattern | cut -d"_" -f2)		#60170515T162341
 dateText=$(echo $datetime | cut -d"T" -f1)		#60170515
 
+#echo "productType " $productType
+echo $tile
+echo $datetime
+echo $dateText
+
 if [ $BYTILE -ne "0" ]; then
+	echo "BY TILE"
 	outDir=$JP2DIR$tile"/"$dateText
 	if [ ! -d $JP2DIR$tile"/" ]; then
 		mkdir $JP2DIR$tile"/"
 	fi
 else
+	echo "BY DATE"
 	outDir=$JP2DIR$dateText
 fi
 

tools/mergeImages.sh

@@ -16,11 +16,8 @@ cd $INDIR
 #BOX=$(polygonToBox.py "$POLYGON")
 BOX=$POLYGON
 ##------------------------------------------------------------------------------
-
-
 echo $DIRNAME
 for f in $(ls);  do echo ${f}; done;
-#exit
 ##------------------------------------------------------------------------------
 MERGEDIMAGE=$DIRNAME"_TCI_10m_merged.tif"
 if [  ! -e $MERGEDIMAGE  ]; then

tools/mergeImages.sh.save

-#!/bin/sh
-DIR=$1	#Output filename
-POLYGON=$2	#Crop Window
-
-echo $DIR
-echo $POLYGON
-
-
-#BOX=$(polygonToBox.py "$POLYGON")
-BOX=$POLYGON
-##------------------------------------------------------------------------------
-cd $DIR
-DIRNAME=$(echo $DIR | cut -d"/" -f1)
-echo $DIRNAME
-##------------------------------------------------------------------------------
-MERGEDIMAGE=$DIRNAME"_TCI_60m_merged.tif"
-if [  ! -e $MERGEDIMAGE  ]; then
-	echo "gdal_merge.py -o "$MERGEDIMAGE $(ls *TCI_10m.jp2)
-	gdal_merge.py -o $MERGEDIMAGE $(ls *TCI_10m.jp2)
-else
-	echo PASSING $MERGEDIMAGE FOUND
-fi
-
-CROPEDIMAGE=../merge_out/$DIRNAME"_TCI_60m.jp2"
-if [  ! -e $CROPEDIMAGE  ]; then
-	gdal_translate -projwin $BOX -projwin_srs WGS84 -ot Byte -of JPEG $MERGEDIMAGE $CROPEDIMAGE
-else
-	echo PASSING $CROPEDIMAGE FOUND
-fi
-##------------------------------------------------------------------------------
-MERGEDIMAGE=$DIRNAME"_SCL_60m_merged.tif"
-if [  ! -e $MERGEDIMAGE  ]; then
-	echo "gdal_merge.py -o "$MERGEDIMAGE $(ls *SCL_60m.jp2)
-	gdal_merge.py -o $MERGEDIMAGE $(ls *SCL_60m.jp2)
-else
-	echo PASSING $MERGEDIMAGE FOUND
-fi
-
-CROPEDIMAGE=../merge_out/$DIRNAME"_SCL_60m.jp2"
-if [  ! -e $CROPEDIMAGE  ]; then
-	gdal_translate -projwin $BOX -projwin_srs WGS84 -ot Byte -of JPEG $MERGEDIMAGE $CROPEDIMAGE
-else
-	echo PASSING $CROPEDIMAGE FOUND
-fi
-
-cd ..
-
-
-

tools/mergeImagesByDirectory.sh

@@ -3,12 +3,15 @@ JP2DIR=$1	#JP2 Directory
 OUTDIR=$2
 POLYGON=$3	#Crop Window
 JOBS=${4:-1}
-DIRNAME=$(echo $JP2DIR | cut -d"/" -f1)
+#DIRNAME=$(echo $JP2DIR | cut -d"/" -f1)
 #echo "POLYGON"
 #echo $POLYGON
+
+
 BOX=$(polygonToBox.py "$POLYGON")
-echo "Polygon BOX"
-echo $BOX
+
+echo "Polygon BOX" $BOX
+
 if [ "$#" -le 1 ]; then
 	echo "Usage " $0 " <Parent Directory with subdirectories> <POLYGON> [jobs]" 
 	exit