Elongated corner reflector?

jmfriedt · August 15, 2022, 9:33am

I am afraid this is trivial but I am unable to find the answer: when mapping corner reflectors using SNAP for interferometric analysis I get elongated shapes rather than the point-like reflectors seen on Sentinel Hub EO Browser
My flowchart is as follows (Makefile calling gpt):

GPT=$(HOME)/snap/bin/gpt
all: target_final.dim
FILENAME1=S1A_IW_SLC__1SDV_20220706T061544_20220706T061612_043976_053FD2_D86C.zip
FILENAME2=S1A_IW_SLC__1SDV_20220718T061545_20220718T061613_044151_054514_02A7.zip

target_split1.dim: $(FILENAME1)
        $(GPT) TOPSAR-Split -Ssource=$(FILENAME1) -PfirstBurstIndex=6 -PlastBurstIndex=6 -Psubswath=IW3 -t target_split1
target_split2.dim: $(FILENAME2)
        $(GPT) TOPSAR-Split -Ssource=$(FILENAME2) -PfirstBurstIndex=6 -PlastBurstIndex=6 -Psubswath=IW3 -t target_split2
target_orb1.dim: target_split1.dim
        $(GPT) Apply-Orbit-File -Ssource=target_split1.dim -t target_orb1
target_orb2.dim: target_split2.dim
        $(GPT) Apply-Orbit-File -Ssource=target_split2.dim -t target_orb2
target_stack.dim: target_orb1.dim target_orb2.dim
        $(GPT) Back-Geocoding -SsourceProducts=target_orb1.dim target_orb2.dim -PdemName="External DEM" -PexternalDEMFile=DEM.tif -PexternalDEMNoDataValue=0.0 -t target_stack.dim  -c 1024M -q 4 # 8192 8
target_interf.dim: target_stack.dim # target_esd.dim single burst => no ESD
        $(GPT) Interferogram -SsourceProduct=target_stack.dim -t target_interf.dim -PsubtractFlatEarthPhase=true -PdemName="External DEM" -PexternalDEMFile="DEM.tif" -PexternalDEMNoDataValue=0.0
target_deburst.dim: target_interf.dim
        $(GPT) TOPSAR-Deburst -SsourceProduct=target_interf.dim -t target_deburst.dim -c 2048M -q 2 
target_topo.dim: target_deburst.dim
        $(GPT) TopoPhaseRemoval -SsourceProduct=target_deburst.dim -PdemName="External DEM" -PexternalDEMFile="DEM.tif" -PexternalDEMNoDataValue=0.0 -t target_topo.dim
target_ml.dim: target_topo.dim
        $(GPT) Multilook -SsourceProduct=target_topo.dim -t target_ml.dim
target_gold.dim: target_ml.dim
        $(GPT) GoldsteinPhaseFiltering -SsourceProduct=target_ml.dim -t target_gold.dim 
target_final.dim: target_gold.dim
        $(eval list1=`ls target_gold.data/q*img | cut -d\/  -f2 | sed 's/\.img//g' | sed 's/q/Phase/g' | tr '\n' ','`)
        $(eval list2=`ls target_gold.data/q*img | cut -d\/  -f2 | sed 's/\.img//g' | sed 's/q/Intensity/g' | tr '\n' ',' | sed 's/,/_db/g'`)
        $(eval list3=`ls target_gold.data/coh*img | cut -d\/  -f2 | sed 's/\.img//g'`)
        $(GPT) Terrain-Correction -SsourceProduct=target_gold.dim -PdemName="External DEM" -PexternalDEMFile="DEM.tif" -PexternalDEMNoDataValue=0.0 -PmapProjection="EPSG:32633" -PsaveSelectedSourceBand=true -PsourceBands="$(list1)$(list2),$(list3)" -t target_final.dim # -PdemName=GETASSE30 -PoutputComplex=true 
        $(eval nom1=`echo $(FILENAME1) | cut -c 20-25`)
        $(eval nom2=`echo $(FILENAME2) | cut -c 20-25`)
        $(eval polar=`echo $(FILENAME2)| cut -c 16`)
        gdal_translate -of GTiff ./target_final.data/Phase*$(polar)$(polar)*.img final_phase$(polar)$(polar)$(nom1)_$(nom2).tif
        gdal_translate -of GTiff ./target_final.data/coh*$(polar)$(polar)*.img final_coh$(polar)$(polar)$(nom1)_$(nom2).tif
        gdal_translate -of GTiff ./target_final.data/Int*$(polar)$(polar)*.img final_intensity$(polar)$(polar)$(nom1)_$(nom2).tif
        gdal_translate -of GTiff ./target_final.data/Int*$(polar)$(polar)*.img final_intensity$(polar)$(polar)$(nom1)_$(nom2).tif
clean:
        rm -rf target*

so Split → Apply precise orbit → Back Geocoding → Interferometry → Deburst → Topo phase removal → Multilook → Goldstein phase filtering → Range Doppler Terrain Correction resulting in the attached thee corner reflector areas:
Screenshot from 2022-08-15 11-18-52
whereas Sentinel Hub EO Browser shows for both dates are “perfectly” overlapping:
2022-08-15-113636_1024x768_scrot
2022-08-15-113618_1024x768_scrot

I have removed ESD since I am working with a single burst and have made as sure as I could that the registration was correct by displaying the RGB image from stacking.
Any obvious step I am missing?

Thanks

mengdahl · August 16, 2022, 8:04am

Did you terrain correct into square ground-range pixels? If not the higher range resolution could be the reason for the apparent elongation.

jmfriedt · August 16, 2022, 5:56pm

I am not sure my reply is appropriate but gdalinfo says of my DEM dataset:

PROJCRS["ETRS89 / UTM zone 33N",
    CONVERSION["UTM zone 33N",
Origin = (420042.500000000000000,8771347.500000000000000)
Pixel Size = (5.000000000000000,-5.000000000000000)
Corner Coordinates:
Upper Left  (  420042.500, 8771347.500) ( 11d14'50.79"E, 78d59'21.37"N)
Lower Left  (  420042.500, 8750672.500) ( 11d18'31.33"E, 78d48'15.88"N)
Upper Right (  445627.500, 8771347.500) ( 12d26'46.57"E, 79d 0' 5.93"N)
Lower Right (  445627.500, 8750672.500) ( 12d29'16.87"E, 78d48'59.69"N)
Center      (  432835.000, 8761010.000) ( 11d52'20.98"E, 78d54'12.82"N)

so I believe by DEM indeed has square pixels when performing the Terrain Correction.

Thanks

jmfriedt · August 18, 2022, 11:48am

Following your comment about the square DEM, I went back to check with SNAP each individual processing step and up to (and including) TopoPhaseRemoval I get a single high-coherence pixel where the corner reflector is located. The elongated shape appears during Multilook processing, which I found awkward since I was running gpt Multilook -SsourceProduct=... -t output without Az and Range arguments so I would have thought they are 1x1.
I hence went back to basics and checked the role of Multilook that I understand as a smoothing filter which is actually not applied in https://step.esa.int/docs/tutorials/S1TBX%20TOPSAR%20Interferometry%20with%20Sentinel-1%20Tutorial_v2.pdf but in http://step.esa.int/docs/tutorials/command_line_inSAR_processing.pdf

Top after Goldstein filtering but no Multilook, bottom with Goldstein after MultiLook

So my question is: what would be the processing advice when targeting a point-like reflector such as a corner reflector. Are MultiLook (and even Goldstein filtering) advised in this circumstance or is Split → Orb → Back Geocoding → Interferogram → Deburst (as shown in the second page of the latter link) sufficient for corner reflector displacement measurement?

Nevertheless after Terrain-Correction I end up with an elongated reflector, with (top) or without (bottom) MultiLook processing, which sounds a bit obvious considering how the non-ML image looks prior to terrain correction.

Would it actually make sense to perform the phase analysis on the output of the Goldstein filtering prior to TerrainCorrection, since I know where the corner reflector is and do not need the geographical coordinates, and avoid spreading energy with this last processing step?

Thanks.

mengdahl · August 18, 2022, 2:11pm

I’m still guessing the terrain-correction step creates the problem. What does the reflector look like in radar-geometry, before the terrain correction step?

Usually people working with corner reflectors want to preserve the resolution, so no multilooking at all or a minimal amount (for example 4x1 in range and azimuth).

jmfriedt · August 18, 2022, 2:28pm

I believe the images prior to terrain correction are the first one I displayed in the previous post (no multilook, with the corner reflector right in the middle as the sharp bright spot in the correlation map) and second one (with multilook, with the corner reflector appearing as an elongated vertical line).
Thanks for clarifying the impact of multilook when working with CR.

jmfriedt · August 18, 2022, 3:49pm

After reading How to import an external DEM into SNAP? - #19 by Parisa I thought the issue might have been that the DEM I am using is in projected coordinates UTM33N but saving as WGS84 and processing again the whole sequence leads to the same elongated shape after Terrain-Correction. Could it be related to the high latitude (which is the reason why I cannot used SRTM)? although I cannot imagine a relationship between latitude and DEM deformation.

mengdahl · August 18, 2022, 7:15pm

Well, corner-reflectors are not elongated and EO Browser is using SNAP for processing S-1 data, so the problem must be somewhere in your processing. To me it looks like you are erroneously terrain correcting non-square radar geometry pixels into square map projected pixels and the ~5 time higher resolution in slant-range is making those reflectors look elongated.

BTW I’m not sure you should want to use goldstein filter either, the corner-reflector phase is unrelated to the phases of the nearby pixels so I don’t see how any “smoothing” could help.

jmfriedt · August 19, 2022, 5:10am

Thank you for confirming that Goldstein is also not needed, that confirms my new understanding of phase analysis of the point-like target.
I did try running the whole processing with SNAP using the GETASSE30 since working at 79 deg N prevents using SRTM and I obtain the same elongated shape, which now that you state it sounds obvious since the target is square in SAR geometry. I wonder though how Sentinel Hub EO Browser ends up with the fine square shape in geographic coordinates … in the mean time I’ll complete the analysis in SAR geometry where I can find the target and identify the phase in the non-projected framework.
Thank you for your support.

ABraun · August 19, 2022, 7:47am

What happens wheny ou terrain correct this (as suggested without Multi Looking or Filtering) but with Nearest Neighbor resampling of the image?

jmfriedt · August 19, 2022, 8:08am

Here is the result of deburst, i.e. no filtering after Interferometry → Deburst (top), and bottom is the Terrain Correction executed in SNAP with my DEM by selecting Nearest Neighbor on both DEM Resampling and Image Resampling instead of the Bilinear Interpolation. Same result though.

Thanks for the support.

mengdahl · August 19, 2022, 8:32am

Hmm, and what if you skip interferometry and terrain-correct the co-registered bursts?

jmfriedt · August 19, 2022, 8:55am

Very interesting indeed: Orbit → Deburst → Terrain Correction using my DEM with projection to WGS84/UTM33N does yield a single pixel! I highlighted the corner reflector manually (red circle) to make identification easier.

Thanks.

I had checked the registration of the two images and the result looked good to me, but maybe I was mistaken (R, G=master, B=slave):

mengdahl · August 19, 2022, 9:04am

Now you can start adding steps one by one to see where things go wrong.

jmfriedt · August 19, 2022, 11:36am

Got it … the default configuration of SNAP Interferogram formation is to compute the correlation on a 10 range-pixels and 2 azimuth-pixel window size. whereas the default setting of gpt is 10x10

Parameter Options:
  -PcohWinAz=<int>                        Size of coherence estimation window in Azimuth direction
                                          Default value is '10'.
  -PcohWinRg=<int>                        Size of coherence estimation window in Range direction
                                          Default value is '10'.

setting -PcohWinAz=2 in gpt Interferogram leads to square pixels after TerrainCorrection.
So that seems to be solved.

Just a last question to make sure I understand all the processing steps: am I correct in understanding that the Topographic Phase correction is “just” a DEM-generated constant phase term that I can also remove for corner reflector double-difference displacement investigation?
Thanks

mengdahl · August 19, 2022, 1:10pm

I think that should be correct.