cuerpos de agua similar a sen2cor

examples/s1ca.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
+import os
+from distutils.dir_util import copy_tree
+
+# 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
+    
+def create_outputdirectory(directory_name):
+    if not os.path.exists(directory_name):
+        os.makedirs(directory_name)
+
+
+# ###############################  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')
+
+
+# crea directorio de salida
+directory=os.path.dirname(imagefullpath)
+os.chdir(directory)
+
+#create_outputdirectory(image_name)
+ca_directory = image_name + '\\IMG_DATA\\ca'
+create_outputdirectory(ca_directory)
+
+#copia datos de antiguo directorio
+copy_tree(imagefullpath,image_name)
+# inicia preprocesamiento para CA
+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 = ca_directory + "\\" + image_name[:-5] + "_Otsu"
+ProductIO.writeProduct(agua_bin, ifg_output, 'GeoTIFF')
+
+