cuerpos de agua con Otsu

examples/ca_geo.py

+# -*- coding: utf-8 -*-
+"""
+Created on Tue Sep 7 10:30:42 2018
+
+@author: A. Salazar,
+"""
+
+# #########################  Importa bibliotecas ##############################  #
+
+import sys
+from snappy import jpy
+from snappy import ProductIO
+from snappy import GPF
+from snappy import HashMap
+import numpy as np
+from skimage import filters
+
+
+
+# Extrae información sobre el producto Sentinel1 GRD solicitado -- sentinel_1_path
+# Regresa la información del producto -- product
+    
+def read_metadata (sentinel_1_path):    
+ 
+    sentinel_1_metadata = "manifest.safe"
+    #s1prd = "%sSAFE/%s" % (sentinel_1_path, sentinel_1_metadata)
+    s1prd = "%s/%s" % (sentinel_1_path, sentinel_1_metadata)
+    reader = ProductIO.getProductReader("SENTINEL-1")
+    product = reader.readProductNodes(s1prd, None)
+
+    return product
+ 
+# Generara los coeficientes de retrodispersión de las bandas VV 
+# del producto solicitado -- input_product  
+# Regresa con calibración radiométrica -- calibrate_product
+
+def radiometric_calibration(input_product):   
+    parameters = HashMap()
+    parameters.put('auxFile', 'Latest Auxiliary File')
+    parameters.put('outputSigmaBand', True)
+    parameters.put('selectedPolarisations', 'VV')
+    calibrate_product = GPF.createProduct('Calibration', parameters, input_product)
+    
+    return calibrate_product
+
+#Realiza el filtrado del speckle al producto previamente calibrado -- calibrate
+#Regresa el producto filtrado -- speckle    
+ 
+def speckle_filtering (calibrate):
+     
+    parameters = HashMap()    
+    parameters.put('filter', 'Lee')
+    parameters.put('filterSizeX', 7)
+    parameters.put('filterSizeY', 7)
+    parameters.put('dampingFactor', 2)
+    parameters.put('edgeThreshold', 5000.0)
+    parameters.put('estimateENL', True)
+    parameters.put('enl', 1.0)    
+    speckle = GPF.createProduct('Speckle-Filter', parameters, calibrate)
+    
+    return speckle
+ 
+#Encuentra el umbral optimo mediante el método de Otsu 
+#Aplica el umbral encontrado al producto de entrada --speckle
+#Este umbral separa lo que es agua de lo que no es agua    
+#Regresa el producto umbralizado  -- agua_otsu
+#Los valores en blaco representan los cuerpos de agua
+    
+def water_detection_otsu(speckle):
+
+    paramToDB = HashMap()
+    
+    paramToDB.put('sourceBandNames', 'Sigma0_' + 'VV')
+    img1 = GPF.createProduct("LinearToFromdB", paramToDB, speckle)
+    img2 = img1.getBand('Sigma0_VV_db')
+    
+    w = img2.getRasterWidth()
+    h = img2.getRasterHeight()
+  
+    img1_data = np.zeros(w * h, np.float32)   
+    img2.readPixels(0, 0, w, h, img1_data)
+    
+    val = filters.threshold_otsu(img1_data) # Otsu method is applied  
+    
+    threshold = val 
+    expression = "Sigma0_VV_db < %s" % threshold    
+    BandDescriptor = jpy.get_type('org.esa.snap.core.gpf.common.BandMathsOp$BandDescriptor')    
+    targetBand = BandDescriptor()
+    targetBand.name = 'agua_Otsu'
+    targetBand.type = 'Float32'
+    targetBand.expression = expression    
+    targetBands = jpy.array('org.esa.snap.core.gpf.common.BandMathsOp$BandDescriptor', 1)
+    targetBands[0] = targetBand    
+    parameters = HashMap()
+    parameters.put('targetBands', targetBands)    
+    #parameters.put('nodataValueAtSea', True)
+    agua_otsu = GPF.createProduct('BandMaths', parameters, img1)   
+
+    return agua_otsu
+    
+
+
+# ###############################  Programa principal   ##################################  #
+
+if len(sys.argv) != 2:
+    print("usage: %s < .SAFE file > " % sys.argv[0])
+    sys.exit(1)
+
+imagefullpath = sys.argv[1]
+
+
+# prepara salida
+s1_pos = imagefullpath.find('S1')
+image_name = imagefullpath[s1_pos:]
+
+image_name = image_name.replace('GRDH','GRDHCA')
+
+sentinel_1_path = imagefullpath  #[:-5] 
+# Lee metadato
+product = read_metadata (sentinel_1_path)
+# Aplica calibración radiométrica
+calibrate = radiometric_calibration(product)
+# Aplica filtrado speckle
+speckle = speckle_filtering(calibrate)
+# Realiza la detección de cuerpos de agua
+agua_bin = water_detection_otsu(speckle)
+# Salva Geotiff con cuerpos de agua
+ifg_output =  image_name[:-5]
+ProductIO.writeProduct(agua_bin, ifg_output, 'GeoTIFF')
+

tools/GRDHProductListToCA.sh

+#!/bin/sh
+#Obtain CA from GRDH Products
+INPUTDIR=$1	#Top directory with a GRDH Subirectory
+
+
+cd $INPUTDIR/GRDH
+
+if [ ! -d ../CA ]; then
+	mkdir ../CA
+fi
+
+for f in $(find . -type f -name '*.zip')
+do
+	cd $INPUTDIR/CA
+	echo $f
+  
+	Z1=$(echo $f | cut -d "_" -f1,2)
+	Z2=$(echo $f | cut -d "_" -f3,4,5,6,7 | cut -d "." -f1)
+        echo $Z1
+	echo $Z2
+
+	ZIPFILE=$Z1"_GRDHCA_"$Z2.zip
+	echo $ZIPFILE
+	if [ ! -e $ZIPFILE ]; then
+		unzip ../GRDH/$f
+		DIRNAME=$(echo $f | cut -d "." -f1,2)
+		echo $DIRNAME
+		python $HOME/Dev/GeoSentinel/examples/ca_geo.py $DIRNAME.SAFE
+		PCA=$(ls -d S1*_GRDHCA_* -t |head -1)
+                echo $PCA
+		PCA=$(echo $PCA | cut -d "." -f1)
+		echo $PCA 
+		zip -r $PCA.zip $PCA.tif
+
+		rm -r *.tif
+                rm -r *.SAFE
+
+	fi
+	cd $INPUTDIR/GRDH
+#	break
+done
+
+