scripts

scripts/GI_L2ASCL_AreaProcessing.sh

+#!/bin/sh
+USERDIR=$1
+BYTILE=${2:-0}
+
+
+#if [ "$#" -nq 1 ]; then
+#	echo "Usage " $0 " <Processs Directory>" 
+
+#	exit
+#fi
+echo "Starting Process..."
+echo $USERDIR"findProducts.json"
+wkt=$(getJSONparameter.py $USERDIR"findProducts.json" wkt)
+RED='\033[0;31m'
+NC='\033[0m' # No Color
+echo "wkt"
+echo $wkt
+
+echo $USERDIR
+cd $USERDIR
+
+##0.- Create Shape File From WKT
+#echo "${RED}Creating Shape file...${NC}"
+#if [ -d "myshape" ]; then
+#	rm -r myshape
+#fi
+
+#wktToShape.py $USERDIR"findProducts.json" "myshape"
+
+
+##1.- Link L2A products 
+#echo ${RED}"Linking Products..."${NC}
+#if [ ! -d "L2A" ]; then
+#	mkdir L2A
+#fi
+#cd L2A
+#	rm *
+#cd ..
+#yes yes | findProducts.py $USERDIR"findProducts.json"
+
+
+##2.- Extract Images
+#echo ${RED}"Extracting JP2 Images..."${NC}
+#if [ ! -d "jp2" ]; then
+#	mkdir jp2
+#fi
+#if [ $BYTILE -ne "0" ]; then
+#        echo "BY TILE"
+#	L2AProductListExtractData.sh $USERDIR"L2A/" $USERDIR"jp2/" 4 1
+
+#else
+#        echo "BY DATE"
+#	L2AProductListExtractData.sh $USERDIR"L2A/" $USERDIR"jp2/" 4 0
+#fi
+
+##3.- Merge Images
+#echo ${RED}"\nMerging Images..."${NC}
+#if [ $BYTILE -ne "0" ]; then
+#	echo "BY TILE"
+#	cd $USERDIR"jp2/"
+#	if [ ! -d "../out" ]; then
+#		mkdir "../out"
+#	fi
+##	ls -d */ | parallel -q --jobs 1 mergeImagesByDirectory.sh $USERDIR"jp2/"{} "$wtk" 4
+#	for i in $(ls -d */) ; do
+#		if [ ! -d "../out/"$i ]; then
+#			mkdir "../out/"$i
+#		fi
+#		mergeImagesByDirectory.sh $USERDIR"jp2/"$i ../../../out/$i "$wkt"
+
+#	done
+#else
+#	echo "BY DATE"
+#        cd $USERDIR"jp2/"
+#        if [ ! -d "../out" ]; then
+#                mkdir "../out"
+#        fi
+#        mergeImagesByDirectory.sh $USERDIR"jp2/" $USERDIR"out/" "$wkt" 4
+
+#fi
+
+##4.- Mask raster with shape file
+#echo ${RED}"\nApplying Mask..."${NC}
+
+#if [ $BYTILE -ne "0" ]; then
+#        echo "BY TILE"
+#        cd $USERDIR"out/"
+#        for d in $(ls -d */) ; do
+#		cd $d
+#	        if [ ! -d "mask" ]; then
+#        	        mkdir mask
+#        	fi
+
+#	        for i in $(ls *.tif) ; do
+#        	        if [  ! -e mask/$i ]; then
+#                	        gdalwarp -cutline  $USERDIR"myshape/wkt.shp" -crop_to_cutline -dstalpha $i mask/$i
+#	                else
+#        	                echo "Passing " $i
+#                	fi
+#	        done
+#		cd ..
+#        done
+
+#else
+#        echo "BY DATE"
+#	cd $USERDIR"out/"
+#	if [ ! -d "mask" ]; then
+#        	mkdir mask
+#	fi
+
+#        for i in $(ls *.tif) ; do
+#		if [  ! -e mask/$i ]; then
+#			gdalwarp -cutline ../myshape/wkt.shp -crop_to_cutline -dstalpha $i mask/$i
+#		else
+#			echo "Passing " $i
+#		fi
+#	done
+#fi
+
+##5.- Delete images with few data
+#echo ${RED}"\nDeleting Images..."${NC}
+
+#if [ $BYTILE -ne "0" ]; then
+#        echo "BY TILE"
+#        cd $USERDIR"out/"
+#	for i in $(ls -d */) ; do
+#		cd $i"mask"
+#			echo $i
+#	        if [ ! -d "SCL" ]; then
+#        	        mkdir SCL
+#        	fi
+#        	cp *SCL_60m.tif SCL
+#        	cd SCL
+#	       	ls *.tif | parallel --jobs 4 imageMissingData.py {} $USERDIR"findproducts.json" 70
+#	        if [ ! -d "../TCI" ]; then
+#        	        mkdir ../TCI
+#        	fi
+#        	for scl in $(ls *.tif)
+#        	do
+#                	fileprefix=$(echo $scl | (cut -d"_" -f1))
+#                	cp ../$fileprefix"_TCI_10m.tif" ../TCI/$fileprefix"_TCI_10m.tif"
+#        	done
+#	        cd $USERDIR"out/"
+
+#	done
+#else
+#        echo "BY DATE"
+#        cd $USERDIR"out/mask/"
+#	if [ ! -d "SCL" ]; then
+#        	mkdir SCL
+#	fi
+#	cp *SCL_60m.tif SCL
+#	cd SCL
+#	ls *.tif | parallel --jobs 4 imageMissingData.py {} $USERDIR"findProducts.json" 70
+
+#	if [ ! -d "../TCI" ]; then
+#        	mkdir ../TCI
+#	fi
+#	for scl in $(ls *.tif)
+#	do
+#        	fileprefix=$(echo $scl | (cut -d"_" -f1))
+#	        cp ../$fileprefix"_TCI_10m.tif" ../TCI/$fileprefix"_TCI_10m.tif"
+#                cd ../TCI
+#	done
+#fi
+##6.-Creating Thumbnails
+##echo ${RED}"Creating Thumbnails"${NC}
+##if [ $BYTILE -ne "0" ]; then
+##        echo "BY TILE"
+##        cd $USERDIR"out/"
+##	for i in $(ls -d */) ; do
+##                cd $i"mask"
+##		if [ ! -d "SCL/thumbnails" ]; then
+##			mkdir SCL/thumbnails
+##	        fi
+##        	if [ ! -d "TCI/thumbnails" ]; then
+##                	mkdir TCI/thumbnails
+##	        fi
+##	        cd TCI
+##        	mogrify -format jpg -quality 100 *.tif
+##	        mogrify -path ./thumbnails -resize 800x600!  *.jpg
+##		cd thumbnails
+##		ffmpeg -pattern_type glob -i "*.jpg" -r 1 -vcodec libx264 $USERDIR"video.mp4"
+
+##		cd $USERDIR"out/"
+##	done
+##else
+##        echo "BY DATE"
+##        cd $USERDIR"out/mask/"
+##        if [ ! -d "SCL/thumbnails" ]; then
+##                mkdir SCL/thumbnails
+##        fi
+##        if [ ! -d "TCI/thumbnails" ]; then
+##                mkdir TCI/thumbnails
+##        fi
+
+##	cd TCI
+##        mogrify -format jpg -quality 100 *.tif
+##	mogrify -path ./thumbnails -resize 10%  *.jpg
+
+##fi
+
+
+##7.-Extract SCL information
+#echo ${RED}"Extracting SCL information..."${NC}
+
+#if [ $BYTILE -ne "0" ]; then
+#        echo "BY TILE"
+#	cd $USERDIR"out/"
+
+#        for i in $(ls -d */) ; do
+#	        cd ../out/$i"mask/SCL"
+#	        if [ ! -d "json" ]; then
+#        	        mkdir json
+#	        fi
+#        	SCLimageListToJSON.sh $USERDIR"out/"$i"mask/SCL/" $USERDIR"findProducts.json" 6
+#	        mv *.json json
+#        	cd json
+#		tile=$(echo $i | cut -d "/" -f1)
+#	        mergeL2ASCL_JSON.py $USERDIR"out/"$i"mask/SCL/json/" > $USERDIR$tile"_sclData.json"
+#		echo $tile
+#		echo "_sclData.json"
+#		cd $USERDIR"out/"
+#	done
+#else
+#        echo "BY DATE"
+#	cd $USERDIR"out/mask/SCL/"
+#	if [ ! -d "json" ]; then
+#        	mkdir json
+#	fi
+#	GI_SCLimageListToJSON.sh $USERDIR"out/mask/SCL/" $USERDIR"findProducts.json" 6
+#	mv *.json json
+#	cd json
+#	mergeL2ASCL_JSON.py $USERDIR"out/mask/SCL/json/" > $USERDIR"scl_data.json"
+#fi
+#8.-Extract SCL morphology
+echo ${RED}"Extracting SCL morphology..."${NC}
+
+if [ $BYTILE -ne "0" ]; then
+        echo "BY TILE"
+
+else
+        echo "BY DATE"
+	cd $USERDIR"out/mask/SCL/"
+	if [ ! -d "morph" ]; then
+        	mkdir morph
+	fi
+	GI_SCLmorphImageListToJSON.sh $USERDIR"out/mask/SCL/" $USERDIR"findProducts.json" 6
+	mv *.json morph
+	cd morph
+	mergeL2ASCL_JSON.py $USERDIR"out/mask/SCL/morph/" > $USERDIR"scl_morph.json"
+fi
+exit
+##8.- Split SCL
+#echo ${RED}"Spliting SCL image..."${NC}
+#if [ $BYTILE -ne "0" ]; then
+#        echo "BY TILE"
+
+#else
+#        echo "BY DATE"
+#        cd $USERDIR"out/mask/SCL/"
+#        if [ ! -d "split" ]; then
+#		mkdir split
+#	fi
+
+#	
+#        for f in $(ls  *.tif) ; do
+#		echo $f
+#		splitSCL.py $f split/
+#	done
+#fi
+
+#9.-Create Tiles
+#echo ${RED}"Creating Tiles..."${NC}
+#if [ $BYTILE -ne "0" ]; then
+#        echo "BY TILE"
+#        cd $USERDIR"out/"
+#        for d in $(ls -d */) ; do
+#		cd $d"mask/TCI/"
+
+#        	if [ ! -d "tiles" ]; then
+#                	mkdir tiles
+#	        fi
+
+#	        for i in $(ls  *.tif) ;
+#        	do
+#                	cd tiles
+#	                gdal2tiles.py ../$i
+#        	        cd ..
+#	        done
+#        	cd $USERDIR"out/"
+#		cd $d"mask/SCL/"
+#                if [ ! -d "tiles" ]; then
+#                        mkdir tiles
+#                fi
+
+#                for i in $(ls  *.tif) ;
+#                do
+#                        cd tiles
+#                        gdal2tiles.py ../$i
+#                        cd ..
+#                done
+#                cd $USERDIR"out/"
+
+#        done
+#else
+#        echo "BY DATE"
+#        cd $USERDIR"out/mask/TCI/"
+#        if [ ! -d "tiles" ]; then
+#                mkdir tiles
+#        fi
+#	for i in $(ls  *.tif) ;
+#	do
+#		cd tiles
+#		gdal2tiles.py --zoom="10-15" ../$i
+#		cd ..
+#	done
+
+#        cd $USERDIR"out/mask/SCL/split/"
+#        if [ ! -d "tiles" ]; then
+#                mkdir tiles
+#        fi
+#        for i in $(ls  *.tif) ;
+#        do
+#                cd tiles
+#                gdal2tiles.py --zoom="10-15" ../$i
+#                cd ..
+#        done
+
+#	
+#fi
+
+exit
+##9.-Create ZIP file
+#echo "Compresing Files in " $USERDIR "..."
+#cd $USERDIR
+#zip -r myzip . --exclude *.json --exclude *.zip --exclude *thumbnails* --exclude *.tif -x ./jp2\* ./L2A\* \*\*/\*/SCL/\*^C
+#echo "END"
+
+#mv *SCL_60m.tif json
+#cd json

scripts/GI_L2ASCLtoJSON.py

+#!/usr/bin/python
+# -*- coding: utf-8 -*-
+import osgeo.ogr as ogr
+import osgeo.osr as osr
+import numpy as np
+import gdal
+import sys
+import os
+import json
+from collections import Counter
+#from collections import OrderedDict
+
+########### lee archivo de configuración ################
+dirname = os.path.dirname(__file__)
+configfile = os.path.join(dirname, '../config/config.json')
+
+with open(configfile, 'r') as f:
+    config = json.load(f)
+
+
+#print(config['PATHS']['PATH_GEOSENTINEL'])
+
+SENTINEL_PATH = config['PATHS']['PATH_GEOSENTINEL']
+###########################################################
+
+sys.path.append(SENTINEL_PATH)
+
+from geosentinel import rasterWkt
+CLASS_MAP = {0:0, 1:1, 2:2, 3:3, 4:4, 5:5, 6:6, 7:7, 8:8, 9:8, 10:8, 11:9}
+#-------------------------------------------------------------------------------
+def L2ASCLtoDict(filename, wkt, wktsub):
+	''' 
+	L2ASCLtoDict
+	'''
+	classMap = lambda x: [ CLASS_MAP[xx] for xx in x ]
+
+	inputImage = gdal.Open(filename)
+	rows, cols, geotransform = inputImage.RasterYSize, inputImage.RasterXSize, inputImage.GetGeoTransform()
+	data = inputImage.GetRasterBand(1).ReadAsArray(0,0,cols,rows)
+
+	data1 = np.bitwise_and(data, rasterWkt.getPolygonArray(inputImage, wkt))
+
+
+	count1 = Counter([])
+	for r in data:
+		count1 = count1+ Counter(classMap(r))
+	sclDict1 = {str(k):v for k,v in dict(count1).items()}
+
+	if wktsub!="":
+		data2 = np.bitwise_and(data1, np.invert(rasterWkt.getPolygonArray(inputImage, wktsub)))
+
+		count2 = Counter([])
+		for r in data2:
+			count2 = count2 + Counter(classMap(r))
+		sclDict2 = {str(k):v for k,v in dict(count2).items()}
+	else:
+		sclDict2=sclDict1
+	return {"areaTotal":sclDict1, "areaNoConservacion":sclDict2}
+#-------------------------------------------------------------------------------
+def main(argv):
+	if len(sys.argv) != 3:
+		print("Usage: " + argv[0] + "<JSON SCL> <JSON Cfg>")
+	else:
+		jsonFile=open(argv[2]).read()
+		cfg = json.loads(jsonFile)
+		dataDict = L2ASCLtoDict(argv[1], cfg['wkt'], cfg["wkt2"])
+		print(json.dumps(dataDict, sort_keys=True))
+#-------------------------------------------------------------------------------
+if __name__ == "__main__":
+	main(sys.argv)

scripts/GI_SCLimageListToJSON.sh

+#!/bin/sh
+DIR=$1
+CGF=$2
+JOBS=${3:-1}
+
+cd $DIR
+ls *.tif | parallel --jobs $JOBS 'GI_L2ASCLtoJSON.py '{} $2 ' > $(echo '{}' | cut -d"." -f1).json'
+
+

scripts/GI_SCLmorphImageListToJSON.sh

+#!/bin/sh
+DIR=$1
+CGF=$2
+JOBS=${3:-1}
+
+cd $DIR
+ls *.tif | parallel --jobs $JOBS 'GI_getComponentsMorphology.py '{} $2 ' > $(echo '{}' | cut -d"." -f1).json'
+
+

scripts/GI_createDirectories.py

+#!/usr/bin/python3
+# -*- coding: utf-8 -*-
+import sys
+#from myModule import myModule
+import json
+import csv
+import os
+import datetime
+#===============================================================================
+def createDirectories(datafile, cfgfile, level1="CVE_ENT", level2="CVE_MUN"):
+	with open(cfgfile) as cfg_file:
+		cfg = json.load(cfg_file)
+	with open(datafile) as csvfile:
+		data = csv.DictReader(csvfile, delimiter=',')
+		for row in data:
+
+			mun = row[level2]
+			ent = row[level1]
+			dirname = ent.zfill(2)+"-"+mun.zfill(3)
+			print(ent)
+#			for k in cfg.keys():
+#				p[k]=cfg[k]
+			print(os.getcwd())
+			cfg["projectDir"]=os.getcwd()+"/"+dirname+"/"
+#			cfg["endDate"]=datetime.datetime.now().strftime("%Y%m%d")
+			cfg["wkt"]=row["WKT"]
+			cfg["wkt2"]=row["WKT2"]
+			if not os.path.exists(dirname):
+				os.mkdir(dirname)
+			with open(dirname+'/findProducts.json', 'w') as outfile:
+    				json.dump(cfg, outfile)
+#===============================================================================
+def main(argv):
+
+        if len(sys.argv) != 3 and len(sys.argv) != 5:
+                print(len(sys.argv))
+                print ("Usage: "+argv[0] + " <csv> <cfg>")
+        elif(len(sys.argv) == 3 ):
+                createDirectories(argv[1], argv[2])
+        else:
+                createDirectories(argv[1], argv[2], argv[3], argv[4])
+
+
+
+if __name__ == "__main__":
+	main(sys.argv)

scripts/GI_getComponentsMorphology.py

+#!/usr/bin/python3
+# -*- coding: utf-8 -*-
+import sys
+import os
+import json
+import csv
+import gdal
+import numpy as np
+import osgeo.ogr as ogr, osr
+sys.path.append("/home/geoint/GeoSentinel")
+from geosentinel import rasterWkt
+from skimage import measure
+AREAm = 60*60 #m²
+AREAkm = AREAm/1000000#km²
+DISTANCE = 60#m
+
+#===============================================================================
+def getComponentsMorphology(image_filename, wkt, wktsub):
+	
+	inputImage = gdal.Open(image_filename)
+	rows, cols, geotransform = inputImage.RasterYSize, inputImage.RasterXSize, inputImage.GetGeoTransform()
+	data = inputImage.GetRasterBand(1).ReadAsArray(0,0,cols,rows)
+
+	data1 = np.bitwise_and(data, rasterWkt.getPolygonArray(inputImage, wkt))
+	if wktsub!="":
+		data2 = np.bitwise_and(data1, np.invert(rasterWkt.getPolygonArray(inputImage, wktsub)))
+	else:
+		data2 = data1
+
+	geometry = ogr.CreateGeometryFromWkt(wkt)
+	inSpatialRef = osr.SpatialReference()
+	inSpatialRef.ImportFromEPSG(4326)
+	outSpatialRef = osr.SpatialReference()
+	outSpatialRef.ImportFromEPSG(32514)
+
+	coordTrans = osr.CoordinateTransformation(inSpatialRef, outSpatialRef)
+	geometry.Transform(coordTrans)
+	wktArea = geometry.GetArea()
+	classDict={}
+	for i in range(12):
+
+		myDict = dict()		
+		dataC = data2==i
+		myDict["area_tot"]=wktArea/1000000
+		myDict["area_verde_tot"] = np.count_nonzero(dataC)
+
+		if myDict["area_verde_tot"]==0.0:
+			continue
+		myDict["area_verde_tot"]=myDict["area_verde_tot"]*AREAkm
+		components_labels = measure.label(dataC)
+		myDict["nAreas"] = int( np.max(components_labels))
+		properties = measure.regionprops(components_labels)
+		
+
+		myDict["gf1"] = 1 / (myDict["area_verde_tot"]/myDict["nAreas"])
+#		myDict["gf2"]
+		myDict["gdv"]=myDict["area_verde_tot"]/myDict["area_tot"]
+		myDict["vp"]=myDict["area_verde_tot"]
+		myDict["tp"]=myDict["area_verde_tot"]/myDict["nAreas"]
+		myDict["fp"]=sum([ p["perimeter"]/p["area"] for p in properties])/myDict["nAreas"]
+
+		myDict["alp"]=sum([ 0 if p["minor_axis_length"]==0 else p["major_axis_length"]/p["minor_axis_length"] for p in properties])/myDict["nAreas"]
+		classDict[str(i)]=myDict
+#		print(myDict)
+#	print(classDict)
+	return classDict
+#===============================================================================
+def main(argv):
+
+	if len(sys.argv) != 3:
+		print (argv[0] + " <image.tif> <cfg.json>")
+	else:
+		with open(argv[2]) as jsonFile:
+			cfg = json.load(jsonFile)
+			dataDict = getComponentsMorphology(argv[1], cfg["wkt"], cfg["wkt2"])
+			print(json.dumps(dataDict, sort_keys=True))
+
+if __name__ == "__main__":
+	main(sys.argv)

scripts/GI_getComponentsMorphology_csv.py

+#!/usr/bin/python3
+# -*- coding: utf-8 -*-
+import sys
+import os
+import json
+import csv
+import gdal
+import numpy as np
+import osgeo.ogr as ogr, osr
+sys.path.append("/home/mchc/git/GeoSentinel")
+from geosentinel import rasterWkt
+from skimage import measure
+AREAm = 60*60 #m²
+AREAkm = AREAm/1000000#km²
+DISTANCE = 60#m
+
+#===============================================================================
+def getComponentsMorphology(image_filename, wktlist):
+	
+	inputImage = gdal.Open(image_filename)
+	rows, cols, geotransform = inputImage.RasterYSize, inputImage.RasterXSize, inputImage.GetGeoTransform()
+	data = inputImage.GetRasterBand(1).ReadAsArray(0,0,cols,rows)
+
+	areaList = []
+	for key, value in wktlist.items():
+		print(key)
+		dataMask = np.bitwise_and(data, rasterWkt.getPolygonArray(inputImage, value))
+		geometry = ogr.CreateGeometryFromWkt(value)
+
+		inSpatialRef = osr.SpatialReference()
+		inSpatialRef.ImportFromEPSG(4326)
+		outSpatialRef = osr.SpatialReference()
+		outSpatialRef.ImportFromEPSG(32514)
+
+		coordTrans = osr.CoordinateTransformation(inSpatialRef, outSpatialRef)
+
+		geometry.Transform(coordTrans)
+		wktArea = geometry.GetArea()
+
+		print(key, wktArea/1000000)
+
+
+		myDict = dict()
+		myDict["clave"] = key
+		dataC = dataMask!=0#==i
+		myDict["area_tot"]=wktArea/1000000
+		myDict["area_verde_tot"] = np.count_nonzero(dataC)
+
+		if myDict["area_verde_tot"]==0.0:
+			continue
+		myDict["area_verde_tot"]=myDict["area_verde_tot"]*AREAkm
+		components_labels = measure.label(dataC)
+		myDict["nAreas"] = int( np.max(components_labels))
+		properties = measure.regionprops(components_labels)
+
+
+		myDict["gf1"] = 1 / (myDict["area_verde_tot"]/myDict["nAreas"])
+##		myDict["gf2"]
+		myDict["gdv"]=myDict["area_verde_tot"]/myDict["area_tot"]
+		myDict["vp"]=myDict["area_verde_tot"]
+		myDict["tp"]=myDict["area_verde_tot"]/myDict["nAreas"]
+		myDict["fp"]=sum([ p["perimeter"]/p["area"] for p in properties])/myDict["nAreas"]
+##		print("axis", [ p["minor_axis_length"] for p in properties])
+		myDict["alp"]=sum([ 0 if p["minor_axis_length"]==0 else p["major_axis_length"]/p["minor_axis_length"] for p in properties])/myDict["nAreas"]
+#		classDict[str(i)]=myDict
+		areaList.append(myDict)
+		print(myDict)
+	with open('morphology.csv', 'w', newline='') as csvfile:
+		writer = csv.DictWriter(csvfile, fieldnames=areaList[0].keys())
+		writer.writeheader()
+		for i in areaList:
+			writer.writerow(i)
+
+##	print(classDict)
+#	return classDict
+#===============================================================================
+def main(argv):
+
+	if len(sys.argv) != 3:
+		print ("Usage: " + argv[0] + " <reference_image.tif> <wkt_list.csv>")
+	else:
+		with open(argv[2]) as csvfile:
+			reader = csv.DictReader(csvfile)
+			areasDic = { row["CVEGEO"]:row["WKT"] for row in reader}
+			
+			dataDict = getComponentsMorphology(argv[1], areasDic)
+
+if __name__ == "__main__":
+	main(sys.argv)

scripts/GI_getMinMax.py

+#!/usr/bin/python3
+# -*- coding: utf-8 -*-
+import sys
+import os
+import json
+import csv
+#===============================================================================
+def getMinMax(clase):
+
+	polylist = []
+	morphList = []
+	for item in os.listdir("."):
+		if os.path.isdir(item):
+			print(item)
+			with open(os.getcwd()+"/"+item+"/scl_data.json") as json_file:
+				dataSCL = json.load(json_file)
+
+			with open(os.getcwd()+"/"+item+"/scl_morph.json") as json_file:
+				dataMorph = json.load(json_file)
+				
+			with open(os.getcwd()+"/"+item+"/findProducts.json") as json_file:
+				dataWKT = json.load(json_file)
+			wkt = dataWKT["wkt"]
+			
+
+			polyDict = {"clave": item, "wkt":wkt}
+
+			cDict = {k:v["areaTotal"][clase] for k, v in dataSCL.items() }
+			maxkey = max(cDict, key=cDict.get )
+			maxval = dataSCL[maxkey]["areaTotal"][clase]
+			polyDict["areaTotal"] = maxval
+			
+			
+			cDict = {k:v["areaNoConservacion"][clase] for k, v in dataSCL.items() if clase in v["areaNoConservacion"]}
+
+			if len(cDict) >0 and clase in dataSCL[maxkey]["areaNoConservacion"]:
+				maxkey = max(cDict, key=cDict.get )
+				maxval = dataSCL[maxkey]["areaNoConservacion"][clase]
+
+			else:
+				maxval=0
+			polyDict["areaNoConservacion"] = maxval
+			polylist.append(polyDict)
+
+
+			tmpDict = dataMorph[maxkey][clase]
+			tmpDict["wkt"] = wkt
+			tmpDict["clave"]=item
+			morphList.append(dataMorph[maxkey][clase])
+	
+	with open('minmax_'+clase+'.csv', 'w') as outfile:
+		writer = csv.DictWriter(outfile, fieldnames=polylist[0].keys())
+		writer.writeheader()
+		for item in polylist:
+			writer.writerow(item)
+
+	with open('morphology_'+clase+'.csv', 'w') as outfile:
+		writer = csv.DictWriter(outfile, fieldnames=morphList[0].keys())
+		writer.writeheader()
+		for item in morphList:
+			writer.writerow(item)
+
+#===============================================================================
+def main(argv):
+
+	if len(sys.argv) != 2:
+		print ("Usage text")
+	else:
+		getMinMax(argv[1])
+
+
+if __name__ == "__main__":
+	main(sys.argv)

scripts/GI_rasterShp.py

+#!/usr/bin/python
+# -*- coding: utf-8 -*-
+#"POLYGON((-89.65513229370117 21.048938357786966,-89.62852478027342 21.047816903874505,-89.6261215209961 21.03339745836918,-89.6465492248535 21.022822311328852,-89.65873718261717 21.039325621609095,-89.65513229370117 21.048938357786966))"
+import osgeo.ogr as ogr
+import osgeo.osr as osr
+import numpy as np
+import gdal
+import sys
+from collections import Counter
+#-------------------------------------------------------------------------------
+def getPolygonArray(image, shp_filename):
+	rows = image.RasterYSize
+	cols = image.RasterXSize
+	geotransform =  image.GetGeoTransform()
+	projection = image.GetProjection()
+	sr = osr.SpatialReference(wkt=projection)
+
+	sr4326 = osr.SpatialReference()
+	sr4326.ImportFromEPSG(4326)
+
+
+	driverMEM	= gdal.GetDriverByName('MEM')
+	driverMemory	= ogr.GetDriverByName('Memory')
+	driverTiff	= gdal.GetDriverByName('GTiff')
+	driverShape	= ogr.GetDriverByName("ESRI Shapefile")
+
+	#RASTER
+	target_ds =  driverMEM.Create('raster.tif', cols, rows, 1, gdal.GDT_Byte) #MEM
+	target_ds.GetRasterBand(1).SetNoDataValue(0)
+	target_ds.GetRasterBand(1).WriteArray(np.zeros((rows, cols), dtype=np.uint8))
+	target_ds.SetProjection(projection)
+	target_ds.SetGeoTransform(geotransform)
+	target_ds.GetRasterBand(1).FlushCache()
+
+	#VECTOR
+	dataSource = driverShape.Open(shp_filename, 0)
+	# Check to see if shapefile is found.
+	if dataSource is None:
+	    print('Could not open', shp_filename)
+	else:
+	    print('Opened %s', shp_filename)
+	    layer = dataSource.GetLayer()
+	    featureCount = layer.GetFeatureCount()
+	    print("Number of features in", shp_filename, featureCount)
+#	rast_ogr_ds = driverMemory.CreateDataSource( 'myshape.shp' )
+#	rast_mem_lyr = rast_ogr_ds.CreateLayer( 'POLYGON', sr4326,  ogr.wkbPolygon)
+#	feat = ogr.Feature( rast_mem_lyr.GetLayerDefn() )
+#	feat.SetGeometry( ogr.CreateGeometryFromWkt(wkt_geom) )
+#	rast_mem_lyr.CreateFeature( feat )
+
+
+	err = gdal.RasterizeLayer(target_ds, [1], layer, burn_values=[255])
+
+	if err != 0:
+		print(err)
+		gdaltest.post_reason( 'got non-zero result code from RasterizeLayer' )
+		return 'fail'
+
+	target_ds.GetRasterBand(1).FlushCache()
+	np_polygon = target_ds.GetRasterBand(1).ReadAsArray(0,0,cols,rows)
+	count = Counter([])
+#	for r in np_polygon:
+#		count = count + Counter(r)
+#	colorDict = {str(k):v for k,v in dict(count).items()}
+#	print(colorDict)
+	target_ds = None
+
+	return np_polygon
+#-------------------------------------------------------------------------------
+def rasterWkt(shpfile, tifffile, outputfile):
+	''' 
+	Draw WKT Polygon into a Image
+	'''
+
+	inputImage = gdal.Open(tifffile)
+	rows = inputImage.RasterYSize
+	cols = inputImage.RasterXSize
+	geotransform =  inputImage.GetGeoTransform()
+	projection = inputImage.GetProjection()
+
+	sr = osr.SpatialReference()
+	sr.ImportFromWkt(projection)
+
+	driverTiff = gdal.GetDriverByName('GTiff')
+	output = driverTiff.Create(outputfile, cols, rows, inputImage.RasterCount, gdal.GDT_Byte)
+
+	np_polygon = getPolygonArray(inputImage, shpfile)
+
+	for b in range(1, output.RasterCount+1):
+		np_input = inputImage.GetRasterBand(b).ReadAsArray(0,0,cols,rows)
+		output.GetRasterBand(b).SetNoDataValue(0)
+		output.GetRasterBand(b).WriteArray(np_polygon)#np.bitwise_and(np_polygon, np_input))
+		output.GetRasterBand(b).FlushCache()
+
+	output.SetGeoTransform(geotransform)
+	output.SetProjection(projection)
+
+	#Close files
+	polygonRaster = None
+	output = None
+	inputImage = None
+
+def main(argv):
+	rasterWkt(argv[1], argv[2], argv[3])
+
+if __name__ == "__main__":
+	main(sys.argv)

scripts/GI_runAnalisis.sh

+#!/bin/sh
+CSV=$1
+CFG=$2
+DIR=$3
+LEVEL1=${4:-CVE_ENT}
+LEVEL2=${5:-CVE_MUN}
+
+if [ ! -d $DIR ]; then
+        mkdir $DIR
+fi
+cd $DIR
+GI_createDirectories.py ../$CSV ../$CFG $LEVEL1 $LEVEL2
+
+for d in $(ls -d */) ; do
+	echo "PROCESSING " $d ...
+	if [ ! -f $(pwd)/"$d/scl_data.json" ] || [ ! -f $(pwd)/"$d/scl_morph.json" ]; then
+		GI_L2ASCL_AreaProcessing.sh $(pwd)/"$d" 0 
+	else 
+	    echo "PASSING " $d
+	fi 
+done
+

scripts/geoJsonToJsonWKT.py

+#!/usr/bin/python3
+# -*- coding: utf-8 -*-
+import sys
+import json
+#from myModule import myModule
+#===============================================================================
+def convertToWKT(filename, outfile):
+	newdata = []
+	with open(filename) as json_file:
+		data = json.load(json_file)
+		for d in data["features"]:
+			prop = d["properties"]
+			print("POLYGON: ", len(d))
+			i = 0
+			if len(d["geometry"]["coordinates"][0]) > 1:
+				print(i)
+				wkt = "POLYGON (( "
+				for p in d["geometry"]["coordinates"][0]:
+					wkt+=str(p[0])+" "+ str(p[1]) +", "
+				i+=1
+				wkt=wkt[:-2]+" ))"
+				print(wkt)
+				
+				newdata.append({"wkt":wkt, "properties":prop})
+
+	with open(outfile, 'w') as outfile:
+		json.dump(newdata, outfile)
+#===============================================================================
+def main(argv):
+
+	if len(sys.argv) != 3:
+		print ("Usage text")
+	else:
+		convertToWKT(argv[1], argv[2])
+
+
+if __name__ == "__main__":
+	main(sys.argv)