S1A product Terrain correction frustration

Dear all,

I am pretty new in the forum and in the SAR topic and need some help here. I am trying to process the following product:

S1A_EW_GRDM_1SDH_20170222T082410_20170222T082510_015394_01942C_945D.SAFE

My Regions Of Interest (ROI) are mainly polar regions where there is no DEM available. Usually those ROIs contain only ocean or is a mixture of ocean and land. I would like to perform the following process:

1 Radiometric correction
2 Speckle filter
3 Terrain correction
4 Mosaicking

(Please let me know at this point if more process needs to be done for the GRD products)

As far as I understand Terrain correction should do also a re-projection to a targeted CRS? My problem is that I am not sure what I miss to perform a correct Terrain correction or/and re-projection. Here is the code:

import snappy
from snappy import ProductIO
from snappy import HashMap
 
import os
from snappy import GPF
from pyproj import Proj
import osgeo.osr

filename  =  S1A_EW_GRDM_1SDH_20170222T082410_20170222T082510_015394_01942C_945D.SAFE
sentinel_1 = ProductIO.readProduct(filename + "/manifest.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)

polarization = 'HH'

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)  

      calib = output + date + "_calibrate_" + polarization 
      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):
    parameters = HashMap() 
    parameters.put('PdemName', 'ASTER 1Sec GDEM')
    parameters.put('imgResamplingMethod', 'BILINEAR_INTERPOLATION')
    parameters.put('incidenceAngleForSigma0', 'Use incidence angle from Ellipsoid')
    parameters.put('PmapProjection', crs)
    output = GPF.createProduct('Terrain-Correction', parameters, source)
    return output

calibrated = do_calibration(sentinel_1, polarization)
filtered = do_speckle_filtering(calibrated)
tercorrected = do_terrain_correction(filtered, srs.ExportToWkt())

The error i am having is:

RuntimeError: org.esa.snap.core.gpf.OperatorException: Entire image is outside of SRTM valid area.
Please use another DEM.

It seems that I can not change the default SRTM DEM. From what I read the ASTER DEM extends up to polar regions which might can be used for my studies.

Could you please help me either to see what I am doing wrong here or, to provide me with a more appropriate method to deal with that? I aim to perform an image enhancement and then to build a mosaic with SAR images.

Many thanks in advance for your time.

Maybe it is better to split your question. Ask the more general question in the s1tbx category. I don’t thing that the SAR experts all look into the python category.
I can’t help you with the SAR questions but maybe with the DEM question.
I see that that you specify for the Terrain Correction the parameter names ‘PdemName’ and ‘PmapProjection’. Remove the ‘P’ from both. Then they might be considered as valid parameter.

Many thanks! It seems this ‘P’ was the problem. I will split the question and add the first part to the SAR experts. Thanks again!

Dear Marco,

After running the corrected version of the code and saving the product as geotiff file:

ProductIO.writeProduct(tercorrected, 'test', 'GeoTIFF')

I ended up with a file of 1.1G size (input is ~200MB) and contains only nan values. Any idea what am I doing wrong? Probably I missed to declare some parameters?

I’m sorry but I don’t know how the radar operators need to be configured.
I’m more a optical guy.

I believe your parameter names are incorrect, inside snappy you don’t prefix them with ‘p’

Here is an example of the terrain correction I do:

params = HashMap()
params.put('demResamplingMethod', 'BILINEAR_INTERPOLATION')
params.put('imgResamplingMethod', 'BILINEAR_INTERPOLATION')
params.put('demName', 'SRTM 3Sec')
params.put('pixelSpacingInMeter', 10.0)
params.put('sourceBands', 'Sigma0_VH')
params.put('nodataValueAtSea', False)

terrain_corrected_product = GPF.createProduct("Terrain-Correction", params, product)
terrain_corrected_product_file_path = get_terrain_path(product_name) + product_name + "_terrain_correction"

ProductIO.writeProduct(terrain_corrected_product,
                       terrain_corrected_product_file_path,
                       "GeoTIFF-BigTIFF")

Woops didn’t notice the comments about ‘P’ beforehand, maybe the code above might still help.

Dear CiaranEvans,

thank you very much for your feedback. I got confused when I called the gpt -h cause it prints the parameters like that : -PnodataValueAtSea= . But now its clear.

If you have some time I would need some more help to understand some more things.

I guess that the “pixelSpacingInMeter” should be the desired spatial resolution of the pixels? The naming is not really self explanatory here :slight_smile: .

You used here "demResamplingMethod"and “imgResamplingMethod”. Please correct me if I am wrong; we use that because probably the pixel center of our targeted output might not be the same with the pixel center of the input. This should be applied though if and only if we change the projection or? If this is the case, shouldn’t you specify also the projection here and put it in the parameters?

Many thanks for your help and your time.

Hi @pakoun

I agree, the naming of the parameters is quite frustrating when trying to duplicate the processes in SNAP/gpt.

You’re correct about ‘pixelSpacingInMeter’ being the desired resolution, unfortunately I was tasked with duplicating an already complete process by someone else at work who uses SNAP desktop so I’m not really able to answer your questions about the resampling, I will endeavour to find you an answer today however!

As far as projection goes, I don’t think it’s required. Having processed a lot of products in the past few weeks I’ve not seen any issues with projections when I’ve put the GeoTiffs into QGIS, it’s lined up perfectly with any other data I’ve had.

I’ll message my colleague regarding the resampling now.

No worries :slight_smile:

Ciaran

Hi @CiaranEvans,

it is very interesting that you are ending up with a perfect match. I can not figure out what I am doing wrong. If you have some spare time might be interesting to take a very quick look on the issue http://forum.step.esa.int/t/python-mosaicing-s1-sar-products-problem-s/5138
The picture is not zoomed in but still probably you can see the misalignment with the coast (around 6k i measured it). I am wondering if this is the case cause my Area Of Interest is very north and probably DEM are not much of a help. May be also something stupid I am doing in the parameters as well.

Many thanks for your help :slight_smile:

Hey @pakoun

Very strange regarding your misalignment, do you mind saying which product file you’re using? Maybe I can run it through my script to see if I get the same issues!

Maybe the shape file (I’m assuming that’s what your land is in the image) is not as accurate as you need? But I doubt it’d be that badly out.

Let me know which product that is and I’ll process it my end!

That would be great and definitely an interesting test. Then I can be sure whether I am doing something stupid :slight_smile: The products I am processing are:
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’]

Reference code you can find under the above mentioned link where I am presenting in detail the whole procedure. The thing is that I am using python to do everything and not the GUI…
The shape file is the standard GSHHS at full resolution. Thats also the default for basemap. I have used also for Antarctic regions a very high resolution product which I also see mismatches with the coast… I have tried lot of things including reprojecting the Ground Control Points to the desired CRS and then make the mosaic (that was GDAL try) but still the problem occurs. This might be related due to the fact that close to the pole lats become singular but I am not sure if this is the case and hot to solve this issue.

Lets see how the test looks like :slight_smile:

Do you have a link to those downloads? I couldn’t get them on https://scihub.copernicus.eu/dhus/#/home

Coastline shapefiles can be very inaccurate in Arctic regions. In my experience GSHHS can be several km off there.
You may want to try Openstreetmap data (e.g. from http://openstreetmapdata.com/)
Looking at http://apps.sentinel-hub.com/sentinel-playground (which seems to use OSM as basemap), the match between the map and Sentinel-2 images is good in that area.

Here are the links:

https://scihub.copernicus.eu/dhus/odata/v1/Products(‘0e90541f-c310-4de9-addf-88b264260a06’)/$value

https://scihub.copernicus.eu/dhus/odata/v1/Products(‘b8f8f4af-5ddc-440b-b64b-e8d49cac0db8’)/$value

https://scihub.copernicus.eu/dhus/odata/v1/Products(‘1643ae7a-a4c3-4b3b-bcad-475cfbf99c66’)/$value

Hope that helps

Thank you for pointing this out. I had the same horrible feeling about the shapefiles but could not imagine how much they might differ. I plotted the GSHHS against the shape file from openstreetmap as you suggested.The area you see is northern part of Greenland. It is interesting (and odd!) that at some regions both shapefiles agree and at some (mainly the more north you go) they significantly differ.

Sorry for the late reply, I left my script running and realised my resolutions were wrong for your data! Had a 40GB+ file being created :’)

Here are my results comparing with Google Satellite in QGIS, they look pretty spot on.

For reference, the parameters I used:

params = HashMap()
params.put('demResamplingMethod', 'BILINEAR_INTERPOLATION')
params.put('imgResamplingMethod', 'BILINEAR_INTERPOLATION')
params.put('demName', 'ASTER 1Sec GDEM')
params.put('pixelSpacingInMeter', 40.0)
params.put('sourceBands', 'Sigma0_HH')
params.put('nodataValueAtSea', False)

Yes you are right thats a perfect match! Thank you for getting into so much trouble! As you and css said these ambiguities come from the shapefiles… May I ask which files are you using? The suggestion from css seems to work really good. I added also a plot to show the differences. For the sake of completeness and helping other users I have to mention that the use of DEM at least for those latitudes seems to be irrelevant with the shore misalignment.

Many thanks both of you guys :slight_smile:

1 Like

No worries :slight_smile: Glad you’ve got it sorted!

I just tried on one but I used: https://scihub.copernicus.eu/dhus/odata/v1/Products(‘0e90541f-c310-4de9-addf-88b264260a06’)/%24value

Good luck in your further work!