Python Mosaicing S1 SAR products problem(s)

Hi again everyone,

I am trying to make a mosaic from S1A products. I managed to go a bit further with the whole procedure but now I am getting some strange stuff while mosaicing. I will start with the products and the code:

import sys
import snappy
from snappy import ProductIO
from snappy import HashMap

import os
from snappy import GPF
from pyproj import Proj
import osgeo.osr
import osr
import gdal, gdalconst
import glob

GPF.getDefaultInstance().getOperatorSpiRegistry().loadOperatorSpis()
HashMap = snappy.jpy.get_type('java.util.HashMap')

safe_files = [S1A_EW_GRDM_1SDH_20170222T082610_20170222T082710_015394_01942C_36A4.SAFE',
              S1A_EW_GRDM_1SDH_20170222T082510_20170222T082610_015394_01942C_A081.SAFE',
              S1A_EW_GRDM_1SDH_20170222T082410_20170222T082510_015394_01942C_945D.SAFE']

proj4str = "+proj={0} +lat_0={1} +lon_0={2} " \
           "+lat_ts={3} +ellps={4} +datum={5} +units={6} " \
           "+lon_wrap={7}".format("stere",85,47,85,"WGS84","WGS84","m",0)
srsproj4 = Proj(proj4str)
srs = osgeo.osr.SpatialReference()
srs.ImportFromProj4(srsproj4.srs)

def do_calibration(source, polarization):
"""
radiometric calibration
"""
    parameters = HashMap() 
    parameters.put('outputSigmaBand', True) 
    parameters.put('sourceBands', 'Intensity_' + polarization) 
    parameters.put('selectedPolarisations', polarization) 
    parameters.put('outputImageScaleInDb', False)  

    target_0 = GPF.createProduct("Calibration", parameters, source) 
    return target_0

def do_speckle_filtering(source):
"""

"""
    parameters = HashMap() 
    parameters.put('sourceBandNames', 'Sigma0_HH')
    parameters.put('filter', 'Refined Lee')
    output = GPF.createProduct('Speckle-Filter', parameters, source)
    return output

def do_terrain_correction(source, crs):
"""

"""
    outname = source.getDisplayName()
    outname = outname + '_tc'

    parameters = HashMap() 
    parameters.put('demName', 'GETASSE30')   #ASTER 1Sec GDEM
    parameters.put('demResamplingMethod', 'BILINEAR_INTERPOLATION')
    parameters.put('imgResamplingMethod', 'BILINEAR_INTERPOLATION')
    parameters.put('pixelSpacingInMeter', 150.0)
    parameters.put('sourceBands', 'Sigma0_HH')
    parameters.put('incidenceAngleForSigma0', 'Use incidence angle from Ellipsoid')
    parameters.put('mapProjection', crs)
    parameters.put('nodataValueAtSea', False)
    output = GPF.createProduct('Terrain-Correction', parameters, source)

   return output

def remove_thermal_noise(source, polarization):
"""

"""
    parameters = HashMap() 
    parameters.put('selectedPolarisations', polarization)
    output = GPF.createProduct('ThermalNoiseRemoval', parameters, source)
    return output

def remove_GRD_border_noise(source, polarization):
"""
removes the black areas around the image
Mask no-value pixels for GRD product
"""
    parameters = HashMap() 
    parameters.put('selectedPolarisations', polarization)
    output = GPF.createProduct('Remove-GRD-Border-Noise', parameters, source)

    return output

def do_mosaic(source):
"""

"""
    parameters = HashMap() 
    parameters.put('pixelSize', 150.0)
    parameters.put('resamplingMethod', 'BILINEAR_INTERPOLATION')
    output = GPF.createProduct('SAR-Mosaic', parameters, source)

    return output

# main script

for s1file in safe_files:
    print(s1file)
    sentinel_1 = ProductIO.readProduct(s1file + "/manifest.safe")
    polarization = 'HH'
    TNfiltered = remove_thermal_noise(sentinel_1, polarization)
    BNfiltered = remove_GRD_border_noise(TNfiltered, polarization)
    calibrated = do_calibration(BNfiltered, polarization)
    filtered = do_speckle_filtering(calibrated)
    tercorrected = do_terrain_correction(filtered, srs.ExportToWkt())

    ProductIO.writeProduct(tercorrected, sentinel_1.getName()+'_final', 'GeoTIFF')
    sentinel_1.dispose()
    sentinel_1.closeIO()

mosaic_files = glob.glob(os.path.join('/home/pakoun/workspace/image_correction', '*final*'))
products = []
for fl in mosaic_files:
    products.append(ProductIO.readProduct(s1file + "/manifest.safe"))

mosaic = do_mosaic(products)
ProductIO.writeProduct(mosaic, 'mosaic', 'GeoTIFF')

When running the Mosaic i get the following warning:

first_line_time metadata value is null
last_line_time metadata value is null

Nevertheless it seems that my file is created. Now if I load it in qgis the mosaic is the (somehow) colored one. As you see this is not a mosaic at all but just only 1 acquisition. Although it supposed to be also TC it looks like it is not projected. For comparison I am plotting also another acquisition (black and white picture) which I performed exactly the same analysis procedure but no Mosaicing. Interestingly enough although I performed a GRD_border_noise removal, the black background which corresponds for the NAN values is present. Further there is a misalignment with the shore which is around 6km…

There might be lot of questions in only one post but I though it would be of great help (providing also the code) for others that they are interested to follow the same procedure.

I trully hope for your input cause I am struggling to understand what is wrong here…

1 Like

Maybe your problem is not so much related to python but to S1 processing in general. What happens when you try the processing in SNAP desktop? Maybe the question is better suited in the s1tbx category.

Hi Marco and thank you for the hint. Is there a way to move this topic to the s1tbx or should i open new topic again? I am just trying to avoid duplications…

I can move it for you

When I run it in snap desctop I do not see any warning and the mosaic seems to be good (that is why I though might be python related problem… maybe a mistake in parameters or so…). I see the whole mosaic of the 3 products (not like the colored one I saw in the picture above). I can open it in snap but only in BEAM-DIMAP format. If I save it as tif and then try to open it, I get this error (my apologies for the long text but might be of your help):

java.lang.IndexOutOfBoundsException: pos < flushedPos!
at javax.imageio.stream.FileImageInputStream.seek(FileImageInputStream.java:143)
at it.geosolutions.imageioimpl.plugins.tiff.TIFFNullDecompressor.decodeRaw(TIFFNullDecompressor.java:214)
at it.geosolutions.imageio.plugins.tiff.TIFFDecompressor.decodeRaw(TIFFDecompressor.java:2168)
at it.geosolutions.imageio.plugins.tiff.TIFFDecompressor.decode(TIFFDecompressor.java:2625)
at it.geosolutions.imageioimpl.plugins.tiff.TIFFNullDecompressor.decode(TIFFNullDecompressor.java:174)
at it.geosolutions.imageioimpl.plugins.tiff.TIFFImageReader.decodeTile(TIFFImageReader.java:1692)
at it.geosolutions.imageioimpl.plugins.tiff.TIFFImageReader.read(TIFFImageReader.java:1970)
at it.geosolutions.imageioimpl.plugins.tiff.TIFFRenderedImage.read(TIFFRenderedImage.java:300)
at it.geosolutions.imageioimpl.plugins.tiff.TIFFRenderedImage.getData(TIFFRenderedImage.java:284)
at org.esa.snap.dataio.geotiff.GeoTiffProductReader.readRect(GeoTiffProductReader.java:256)
at org.esa.snap.dataio.geotiff.GeoTiffProductReader.readBandRasterDataImpl(GeoTiffProductReader.java:158)
at org.esa.snap.core.dataio.AbstractProductReader.readBandRasterData(AbstractProductReader.java:250)
at org.esa.snap.core.image.BandOpImage.computeProductData(BandOpImage.java:67)
at org.esa.snap.core.image.RasterDataNodeOpImage.computeRect(RasterDataNodeOpImage.java:127)
at javax.media.jai.SourcelessOpImage.computeTile(SourcelessOpImage.java:137)
at com.sun.media.jai.util.SunTileScheduler.scheduleTile(SunTileScheduler.java:904)
at javax.media.jai.OpImage.getTile(OpImage.java:1129)
at org.esa.snap.core.datamodel.StxFactory.accumulateTile(StxFactory.java:339)
at org.esa.snap.core.datamodel.StxFactory.accumulate(StxFactory.java:323)
at org.esa.snap.core.datamodel.StxFactory.accumulate(StxFactory.java:296)
at org.esa.snap.core.datamodel.StxFactory.create(StxFactory.java:200)
at org.esa.snap.core.datamodel.StxFactory.create(StxFactory.java:274)
at org.esa.snap.core.datamodel.RasterDataNode.computeStxImpl(RasterDataNode.java:2427)
at org.esa.snap.core.datamodel.Band.computeStxImpl(Band.java:504)
at org.esa.snap.core.datamodel.RasterDataNode.getStx(RasterDataNode.java:2394)
at org.esa.snap.core.datamodel.RasterDataNode.createDefaultImageInfo(RasterDataNode.java:1802)
at org.esa.snap.core.datamodel.Band.createDefaultImageInfo(Band.java:483)
at org.esa.snap.core.datamodel.RasterDataNode.getImageInfo(RasterDataNode.java:1784)
at org.esa.snap.core.datamodel.RasterDataNode.getImageInfo(RasterDataNode.java:1767)
at org.esa.snap.core.image.ImageManager.prepareImageInfos(ImageManager.java:834)
at org.esa.snap.core.image.ColoredBandImageMultiLevelSource.create(ColoredBandImageMultiLevelSource.java:51)
at org.esa.snap.core.image.ColoredBandImageMultiLevelSource.create(ColoredBandImageMultiLevelSource.java:47)
at org.esa.snap.core.image.ColoredBandImageMultiLevelSource.create(ColoredBandImageMultiLevelSource.java:41)
at org.esa.snap.ui.product.ProductSceneImage.(ProductSceneImage.java:69)
at org.esa.snap.rcp.actions.window.OpenImageViewAction.createProductSceneImage(OpenImageViewAction.java:271)
at org.esa.snap.rcp.actions.window.OpenImageViewAction.access$100(OpenImageViewAction.java:69)
at org.esa.snap.rcp.actions.window.OpenImageViewAction$1.doInBackground(OpenImageViewAction.java:236)
at org.esa.snap.rcp.actions.window.OpenImageViewAction$1.doInBackground(OpenImageViewAction.java:213)
at com.bc.ceres.swing.progress.ProgressMonitorSwingWorker.doInBackground(ProgressMonitorSwingWorker.java:55)
at javax.swing.SwingWorker$1.call(SwingWorker.java:295)
at java.util.concurrent.FutureTask.run(FutureTask.java:266)
at javax.swing.SwingWorker.run(SwingWorker.java:334)
at java.util.concurrent.ThreadPoolExecutor.runWorker(ThreadPoolExecutor.java:1142)
at java.util.concurrent.ThreadPoolExecutor$Worker.run(ThreadPoolExecutor.java:617)
at java.lang.Thread.run(Thread.java:745)
Caused: java.util.concurrent.ExecutionException
at java.util.concurrent.FutureTask.report(FutureTask.java:122)
at java.util.concurrent.FutureTask.get(FutureTask.java:192)
at javax.swing.SwingWorker.get(SwingWorker.java:602)
[catch] at org.esa.snap.rcp.actions.window.OpenImageViewAction$1.done(OpenImageViewAction.java:221)
at javax.swing.SwingWorker$5.run(SwingWorker.java:737)
at javax.swing.SwingWorker$DoSubmitAccumulativeRunnable.run(SwingWorker.java:832)
at sun.swing.AccumulativeRunnable.run(AccumulativeRunnable.java:112)
at javax.swing.SwingWorker$DoSubmitAccumulativeRunnable.actionPerformed(SwingWorker.java:842)
at javax.swing.Timer.fireActionPerformed(Timer.java:313)
at javax.swing.Timer$DoPostEvent.run(Timer.java:245)
at java.awt.event.InvocationEvent.dispatch(InvocationEvent.java:311)
at java.awt.EventQueue.dispatchEventImpl(EventQueue.java:756)
at java.awt.EventQueue.access$500(EventQueue.java:97)
at java.awt.EventQueue$3.run(EventQueue.java:709)
at java.awt.EventQueue$3.run(EventQueue.java:703)
at java.security.AccessController.doPrivileged(Native Method)
at java.security.ProtectionDomain$JavaSecurityAccessImpl.doIntersectionPrivilege(ProtectionDomain.java:76)
at java.awt.EventQueue.dispatchEvent(EventQueue.java:726)
at org.netbeans.core.TimableEventQueue.dispatchEvent(TimableEventQueue.java:159)
at java.awt.EventDispatchThread.pumpOneEventForFilters(EventDispatchThread.java:201)
at java.awt.EventDispatchThread.pumpEventsForFilter(EventDispatchThread.java:116)
at java.awt.EventDispatchThread.pumpEventsForHierarchy(EventDispatchThread.java:105)
at java.awt.EventDispatchThread.pumpEvents(EventDispatchThread.java:101)
at java.awt.EventDispatchThread.pumpEvents(EventDispatchThread.java:93)
at java.awt.EventDispatchThread.run(EventDispatchThread.java:82)

Many thanks!

Some more insights on how far I am… It seemed that one of the numerous problems I have in mosaicing was related to the band selection (which I was not selecting it explicitly). This is why I could not read the geotiff file into qgis. The thing is that even though i did not select the band the mosaic was running and produced 3 times larger file which you could not open it. I modified the code as following:

parameters = HashMap() 
parameters.put('pixelSize', 250.0)
parameters.put('resamplingMethod', 'BILINEAR_INTERPOLATION')
parameters.put('sourceBands', 'Sigma0_HH')
output = GPF.createProduct('SAR-Mosaic', parameters, source)

The output in qgis is the following:

If I look into the parameter options for the SAR-mosaic I do not see any crucial parameter that i missed to define
Parameter Options:
-Paverage= Average the overlapping areas Default value is ‘true’.
-PconvergenceThreshold= Convergence threshold for Relaxed Gauss-Seidel method
Default value is ‘1e-4’.
-Pfeather= Feather amount around source image Default value is ‘0’.
-PgradientDomainMosaic= Gradient Domain Mosaic Default value is ‘false’.
-PmaxIterations= Maximum number of iterations Default value is ‘5000’.
-PnormalizeByMean= Normalize by Mean Default value is ‘true’.
-PpixelSize= Pixel Size (m) Default value is ‘0’.
-PresamplingMethod= The method to be used when resampling the slave grid onto the master grid. Value must be one of ‘NEAREST_NEIGHBOUR’, ‘BILINEAR_INTERPOLATION’, ‘CUBIC_CONVOLUTION’.
Default value is ‘NEAREST_NEIGHBOUR’.
-PsceneHeight= Target height Default value is ‘0’.
-PsceneWidth= Target width Default value is ‘0’.
-PsourceBands=<string,string,string,…> The list of source bands.

If I run the SNAP desctop load the files which I already performed image enhancement and make mosaic which I save as geotiff then the file can not be loaded at all. I am getting really confused here. Is it something I miss in my mosaic function or something wrong while saving the geotiff file?

Hi @pakoun

Just looking at your QGIS Screenshot, have you tried altering the values for white + black?

I know I had an issue where GeoTiffs would appear black but changing the min and max values displayed the correct 'pixel’s so to speak.

For example here is one of my SAR GeoTiffs with the same min/max as yours:

Here’s my min/max values I was using:

And here’s the resulting display in QGIS:

@pakoun have you also tried using 'GeoTIFF-BigTIFF' as your format for saving in Python?

Hi @CiaranEvans,

yeap I edited the black white values range and what I got was black white noise (which I forgot to paste in my previous post) like you watch a film from 50s or something (if you are a bit older you remember maybe those movies :joy: ) The second weird thing is that from this noise that I retrieved, was showing only 1 out of the 3 products I used. Nevertheless you see here a very wide result. So basically that was not a mosaic. I am just thinking if this is caused due to NAN values around the acquisition (thats the black area when you insert ESA S1 product directly to qgis). If so, somehow those should be depricated; something similar to the gdal command :

gdalwarp  "-srcnodata None"   "-dstnodata 0"

dunno though how to do this with snappy.

The GeoTIFF-BigTIFF which I also tried, somehow creates a huge file (1.5 Gb)
Similar with the normal GeoTIFF format, but now the black area is squared:

Something is seriously wrong here :joy:

Hey @pakoun

Oh no! It’s not looking great! Unfortunately I’m at bit of a loss now, not sure what to recommend from here.

Hopefully someone will be able to help, I’d be interested to know what’s causing it!

Yeap I am also lost… Just to be sure that there is no code problem… do you see something wrong with the terrain correction parameters here? :

parameters = HashMap() 
parameters.put('demName', 'GETASSE30')   #ASTER 1Sec GDEM
parameters.put('demResamplingMethod', 'BILINEAR_INTERPOLATION')
parameters.put('imgResamplingMethod', 'BILINEAR_INTERPOLATION')
parameters.put('pixelSpacingInMeter', 150.0)
parameters.put('sourceBands', 'Sigma0_HH') #Sigma0_HH
parameters.put('incidenceAngleForSigma0', 'Use incidence angle from Ellipsoid')
parameters.put('mapProjection', crs)
parameters.put('nodataValueAtSea', False)
output = GPF.createProduct('Terrain-Correction', parameters, source)

The source product is the outcome of the image enhancement (noise removal calibration etc.)
Many thanks again!

@pakoun the only one which screams out at me is pixelSpacingInMeter I’ve always set it at the default, 150 seems quite big? I’ve always got mine at 10.0 for S1 SAR stuff

well I need to reduce the file size, so usually I am working with up to 250 m resolution :slight_smile: Its definitely coarser but still the spatial structure I want to see is resolved. Consider also that if you want to plot a jpeg or something, even if you have 10m resolution, the dpi that you define when saving a plot most likely, will not be able to resolve the nominal spatial resolution.

Ah I see, understandable. Unfortunately I don’t know then, I don’t use mapProjection or incidenceAngleForSigma0 in my pipeline, but I’m not sure if that’s doing anything to affect it this badly.

When i try to make a mosaic, i have a error:

mosaic_pro = GPF.createProduct(‘SAR-Mosaic’, mosaicking, test_mosaic)
RuntimeError: org.esa.snap.core.gpf.OperatorException: Product ‘S1A_IW_GRDH_1SDV_20170106T2236_mosaic_’ has no geo-coding.

Does anyone know how to fix this error, this is my code:

    mosaicking = HashMap()
    mosaicking.put('pixelSize', 10.0)
    mosaicking.put('resamplingMethod', 'BILINEAR_INTERPOLATION')
    mosaicking.put('sourceBands', 'Sigma0_VV')
    mosaic_output = mosaic_path + substring + "_mosaic_"
    mosaic_pro = GPF.createProduct('SAR-Mosaic', mosaicking, test_mosaic)
    ProductIO.writeProduct(mosaic_pro, mosaic_output, 'BEAM-DIMAP')
    mosaikings.append(mosaic_pro)
    mosaic_pro.dispose()

The problem is the test_mosaic. How do you create it?
The error message says that is does not have a geo-coding. Probably it is not well defined.
It is not a normal S1 product. The _mosaic_ name extension identifies this. So probably the problem is one step before the mosaicing.

Thanks marpet

I read product after Terrain Correction. This is full code i made:

import glob
mosaic_files = “D:\SENTINEL\Results\6.TerrainCorrection\*.dim”
files = glob.glob(mosaic_files)
for mosaic in files:
mosaic_path = “D:\SENTINEL\Results\7.Mosaic\”
mosaic_data = ProductIO.readProduct(mosaic)

mosaikings = []
mosaikings.append(mosaic_data)
names = mosaic_data.getName()
substring = names[:-48]
test_mosaic = []  
if substring in mosaikings:
    print substring
else:
    test_mosaic.append(substring)
    
for test_mosaic in mosaikings:
    mosaicking = HashMap()
    mosaicking.put('pixelSize', 10.0)
    mosaicking.put('resamplingMethod', 'BILINEAR_INTERPOLATION')
    mosaicking.put('sourceBands', 'Sigma0_VV')
    mosaic_output = mosaic_path + substring + "_mosaic_"
    mosaic_pro = GPF.createProduct('SAR-Mosaic', mosaicking, test_mosaic)
    ProductIO.writeProduct(mosaic_pro, mosaic_output, 'BEAM-DIMAP')
    mosaikings.append(mosaic_pro)
    mosaic_pro.dispose()
print '7.FINISH MOSAICKING DATA...........................DONE'

It really seems that one of the products “D:\SENTINEL\Results\6.TerrainCorrection*.dim” has no geo-coding.
Maybe something went wrong during the terrain-correction. Have you opened the S1A_IW_GRDH_1SDV_20170106T2236_mosaic_ in SNAP Desktop and checked if it is geo-coded. You can also try to exclude it and see what happens.

Hi everyone, its been some time since my last efforts. I am trying again to build a SAR mosaic. Four Sentinel1 files are preprocessed with thermal and ground border noise removal, radiometric calibration, speckle filtering, terrain correction and saved as geotiffs with a spatial resolution of 150m, using snappy. Their file size at this stage is around 100Mb each. Gdalinfo shows that float32 is being used. Those files will be used as input for the mosaic. So far so good.
I have the following situation:

When using SNAP and saving it as geotiff the mosaic looks fine with a file size of 1.7 Gb (BTW isn’t a bit large??).

Now I am trying to use snappy. First attempt with GeoTIFF:
java.lang.IllegalStateException: File size too big. TIFF file size is limited to [4294967296] bytes!
That is odd, why with snappy the size seems to be more than 4G?

Next attempt with GeoTIFF-Big-TIFF
After 3h the file size is more than 8G and it seems that the script is hanging…
Attempt to save it as NetCDF-CF produces the following error:
RuntimeError: java.lang.IllegalArgumentException: Variable size in bytes 8545341456 may not exceed 4294967292

Does somehow the ProductIO.writeProduct method implicitly assumes float64 or something? If this is the problem, iIs there a way to set this to different format? Is there any work-around to save it with a geotiff format using snappy?