i used snappy to automate user requested sar-processing for seawall detection, initially i used SNAP-desktop to get processed and it worked fine ,
but when i add all same details the snappy gave entire blank image after terrain correction , any one can help ???
Python Script :
# Need to configure Python to use the SNAP-Python (snappy) interface(https://senbox.atlassian.net/wiki/spaces/SNAP/pages/50855941/Configure+Python+to+use+the+SNAP-Python+snappy+interface)
# Read in unzipped Sentinel-1 GRD products (EW and IW modes)
# Parameters to provide: path, outpath, region_code
import datetime
import time
from snappy import GPF
from snappy import ProductIO
from snappy import HashMap
import os, gc
from snappy import jpy
from snappy import File,ProgressMonitor
from snappy import GPF
import geopandas as gpd
from shapely.geometry import MultiPolygon, Point,Polygon
def createProgressMonitor():
PWPM = jpy.get_type('com.bc.ceres.core.PrintWriterProgressMonitor')
JavaSystem = jpy.get_type('java.lang.System')
monitor = PWPM(JavaSystem.out)
return monitor
def do_apply_orbit_file(source):
print('\tApply orbit file...')
parameters = HashMap()
parameters.put('Apply-Orbit-File', True)
output = GPF.createProduct('Apply-Orbit-File', parameters, source)
return output
def do_thermal_noise_removal(source):
print('\tThermal noise removal...')
parameters = HashMap()
parameters.put('removeThermalNoise', True)
output = GPF.createProduct('ThermalNoiseRemoval', parameters, source)
return output
def do_calibration(source, polarization, pols):
print('\tCalibration...')
parameters = HashMap()
parameters.put('outputSigmaBand', True)
if polarization == 'DH':
parameters.put('sourceBands', 'Intensity_HH,Intensity_HV')
elif polarization == 'DV':
parameters.put('sourceBands', 'Intensity_VH,Intensity_VV')
elif polarization == 'SH' or polarization == 'HH':
parameters.put('sourceBands', 'Intensity_HH')
elif polarization == 'SV':
parameters.put('sourceBands', 'Intensity_VV')
else:
print("different polarization!")
parameters.put('selectedPolarisations', pols)
parameters.put('outputImageScaleInDb', False)
output = GPF.createProduct("Calibration", parameters, source)
return output
def do_speckle_filtering(source):
print('\tSpeckle filtering...')
parameters = HashMap()
parameters.put('filter', 'Lee')
parameters.put('filterSizeX', 5)
parameters.put('filterSizeY', 3)
parameters.put('NumberOfLooks', 4)
output = GPF.createProduct('Speckle-Filter', parameters, source)
return output
def do_terrain_correction(source, proj, downsample):
print('\tTerrain correction...')
parameters = HashMap()
parameters.put('demName', 'SRTM 1Sec Grid')
parameters.put('demResamplingMethod', 'BILINEAR_INTERPOLATION')
parameters.put('imgResamplingMethod', 'BILINEAR_INTERPOLATION')
#parameters.put('mapProjection', proj) # comment this line if no need to convert to UTM/WGS84, default is WGS84
parameters.put('saveProjectedLocalIncidenceAngle',False)
parameters.put('saveSelectedSourceBand', True)
parameters.put('pixelSpacingInMeter', 10.0)
parameters.put('nodataValueAtSea', False)
parameters.put('maskOutAreaWithoutElevation',False)
parameters.put('Latitude&Longitude',True)
output = GPF.createProduct('Terrain-Correction', parameters, source)
return output
def do_subset(source, wkt):
print('\tSubsetting...')
parameters = HashMap()
parameters.put('geoRegion', wkt)
output = GPF.createProduct('Subset', parameters, source)
return output
def save(file,output_path,file_format) :
terrain=output_path
WriteOp = jpy.get_type('org.esa.snap.core.gpf.common.WriteOp')
writeOp = WriteOp(file, File(terrain),file_format)
writeOp.writeProduct(ProgressMonitor.NULL)
def point_to_circle(lat,lon,radius):
radius/=100.0
p=Point([(lon,lat)])
g=gpd.GeoSeries([p])
circle=g.buffer(radius)[0]
coords=list(circle.exterior.coords)
wkt=circle.to_wkt()
return wkt,coords
## All Sentinel-1 data sub folders are located within a super folder (make sure the data is already unzipped and each sub folder name ends with '.SAFE'):
path = 'tmp/'
outpath = 'tmp/'
if not os.path.exists(outpath):
os.makedirs(outpath)
## UTM projection parameters
proj = '''PROJCS["UTM Zone 4 / World Geodetic System 1984",GEOGCS["World Geodetic System 1984",DATUM["World Geodetic System 1984",SPHEROID["WGS 84", 6378137.0, 298.257223563, AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich", 0.0, AUTHORITY["EPSG","8901"]],UNIT["degree", 0.017453292519943295],AXIS["Geodetic longitude", EAST],AXIS["Geodetic latitude", NORTH]],PROJECTION["Transverse_Mercator"],PARAMETER["central_meridian", -159.0],PARAMETER["latitude_of_origin", 0.0],PARAMETER["scale_factor", 0.9996],PARAMETER["false_easting", 500000.0],PARAMETER["false_northing", 0.0],UNIT["m", 1.0],AXIS["Easting", EAST],AXIS["Northing", NORTH]]'''
for folder in os.listdir(path):
gc.enable()
gc.collect()
sentinel_1 = ProductIO.readProduct(path + "/" + folder + "/manifest.safe")
print(sentinel_1)
loopstarttime=str(datetime.datetime.now())
print('Start time:', loopstarttime)
start_time = time.time()
## Extract mode, product type, and polarizations from filename
modestamp = folder.split("_")[1]
productstamp = folder.split("_")[2]
polstamp = folder.split("_")[3]
polarization = polstamp[2:4]
if polarization == 'DV':
pols = 'VH,VV'
elif polarization == 'DH':
pols = 'HH,HV'
elif polarization == 'SH' or polarization == 'HH':
pols = 'HH'
elif polarization == 'SV':
pols = 'VV'
else:
print("Polarization error!")
## Start preprocessing:
#applyorbit = do_apply_orbit_file(sentinel_1)
#wkt='POLYGON ((106.501 20.301, 106.501 21.801,108.001 21.801,108.001 20.301,106.501 20.301))'
wkt,coords=point_to_circle(21.626793, 108.453163,50)
# save(thermaremoved,outpath + '/' + folder[:-5],'BEAM-DIMAP')
# del thermaremoved
# sentinel_1.dispose()
# sentinel_1.closeIO()
# thermaremoved=ProductIO.readProduct(outpath + '/' + folder[:-5]+'.dim')
calibrated = do_calibration(sentinel_1, polarization, pols)
#save(calibrated,outpath + '/' + folder[:-5],'GeoTIFF')
# del calibrated
# thermaremoved.dispose()
# thermaremoved.closeIO()
# calibrated=ProductIO.readProduct(outpath + '/' + folder[:-5]+'.dim')
down_filtered = do_speckle_filtering(calibrated)
# save(down_filtered,outpath + '/' + folder[:-5],'BEAM-DIMAP')
# calibrated.dispose()
# calibrated.closeIO()
# down_filtered=ProductIO.readProduct(outpath + '/' + folder[:-5]+'.dim')
tercorrected = do_terrain_correction(down_filtered, proj, 0)
tercorrected=do_subset(tercorrected, wkt)
#save(tercorrected,'tmp/sar','GeoTIFF')
# down_filtered.dispose()
# down_filtered.closeIO()
# del down_filtered
ProductIO.writeProduct(tercorrected, 'tmp/sar', 'GeoTIFF-BigTIFF')
sentinel_1.dispose()
sentinel_1.closeIO()
print("--- %s seconds ---" % (time.time() - start_time))