Sentinel Interferometric Wide Focusing
LucaMancaITA opened this issue · 18 comments
Dear Jimfriedt,
first of all I would like to thank you again for the amazing work you provided.
After succesfully decoded a raw Sentinel IW image, I am trying to focus it. For what concern the range compression I am using the formula available in your code (in accordance with the Level 1 Detailed Algorithm Definition document) and the result seems to be something similar to what one can expect after a range compression (see figure below).
However I still have issues in the azimuth compression process. My question is: why you did not used the following azimuth chirp formula used in the Range Doppler Algorithm?
C(s) = exp( j * (4* pi / lambda) * R(s)), approximating R(s) using a Taylor expansion
Below the images of the decoded raw IW Sentinel image (swath 2), the image after range compression (still not understand why it is mirrored along the y axis) and the corresponding slc image (read using snappy toolbox)
My other question is: is there a way to focus raw Sentinel SM/IW images without creating the reference azimuth chirp doing a crude search for the strongest target in the scene?
Thanks again for your time,
Luca
Dear Luca Manca,
Did you use the same read_bin.m file as provided in the Code to produce the range compressed image or did some changes to it, Ive been trying to process a Decoded Raw Data file through the provided read_bin.m file but havent been able to do so. Could you kindly guide me?
Thank you,
Best Regards
Hi,
yes I've used exactly the code read_bin.m to perform the range compression step. Could you explain the error that you faced during the process?
Best regards,
Luca
The following is the image plot of the Binary file i-e resultSW11_T1286748902_NQ12061.bin (unfocused as i understand)
I use the following code to plot 'sol(:,k)' as image after range compression to check for some range compression data
% range compression completed
figure (2)
subplot (211)
imagesc(abs(sol(:,k)))
subplot (212)
plot (sol(:,k))
clear t
fclose(f)
I get the following results.
I am not sure if it is being compressed properly, because the display of my results vary alot from yours or jmfried`s. (Please guide if there is any other way)
i get the following result after imagesc function is applied to ssol(:,k) at the end of the m file,
This result does not explain anything as far as i understand.
in the end i apply the imagesc function to the X-matrix obtained after running the Read_bin.m file.
it gives me the following results,
The issue I have is, I believe i am not getting the desired results as can be seen from the plots. the data is of a region over Manaus, near Brazil which comprises of both Rivers and Land.
Can you please guide me towards a better way to plot these or assess the issue, Kindly.
Thank you for your reply,
Best Regards
- LucaMancaITA: your range compression seems to operate nicely. Since your range axis is the X axis, and usually correlation is iFFT(FFT(signal).*FFT(chirp)), axis flipping could be due to an FFT instead of iFFT when returning to space domain after the product of the FFT. Replacing FFT with iFFT should swap the X-axis.
- LucaMancaITA: azimuth compression using PRI and orbital parameters can be done and is described in the literature, but is beyond my understanding (at the moment). I refer to https://topex.ucsd.edu/gmtsar/tar/GMTSAR_2ND_TEX.pdf pp.51--59 when addressing the theoretical approach to azimuth compression, but I do not understand the underlying concepts (yet). In http://jmfriedt.free.fr/glmf_level0.pdf (in French, maybe Google Translate can help) I compare in Fig 17 (page 15) the parabolic phase evolution as obtained from the strongest reflector with the application of the model mentioned above and the match seems quite good to me. But I cannot help to understand where the mathematics come from in this model.
- wiggedcrawdad951: your imagesc are too saturated to see anything. I am not sure I understand what you mean with imagesc(abs(sol(:,k))) (you imagesc a vector which does not seem to make much sense to me) and plot(sol) will indeed plot a constellation of the complex values which also does not seem to make much sense in this context (I understand constellations for digital communication after carrier and symbol synchronization, but not in the context of RADAR processing).
Best wishes and thank you for your interest, JM
Jmfriedt: Right, can you guide me as to how can i get good results, and display range compressed data and the final image properly, in the range compression phase i was trying to plot the range compressed data which i failed to do so as can be seen.
Thank you for your time.
best regards.
Please provide the reference to the exact RAW dataset you are investigating (the filename on the Copernicus repository) so I can download and reproduce your processing.
Thanks.
This is the data set that I have been using,
S1B_IW_RAW__0SDV_20201014T221423_20201014T221455_023814_02D411_C1D3.SAFE
if it is offline in the copernicus respository, it can be accessed from here,
(Alaska satellite Facility)
Thank you, Best Regards
I processed the same dataset and indeed abs(x(1:10:end,:)) after x=read_bin('resultSW11_T1286748902_NQ12061.bin',2000,12061*2); does not show much feature.
However I had a look at the (more recent) S1B_IW_SLC__1SDV_20210729T221428_20210729T221456_028014_03578E_A6B0
over the same area and unless you select a dataset on the southwestern lake, you might not be able to see much feature on the homogeneous Amazon forest. Sao Paulo has been selected as a reference dataset because it features StripMap (SM) measurements (as opposed to Interferometric Wide (IW) -- easier for beginning) and point-like ships waiting in the harbour (actually I have been told that the harbour is not called Sao Paulo but Santos), so unless you have a specific reason for analyzing the raw data of a forest covered area, I would consider at first areas with easier to identify features.
Maybe you will find it more satisfactory to run
x=read_bin('resultSW10_T1286748884_NQ11924.bin',1400,11924*2);
imagesc((abs(x(1:10:end,:))),[0 1000])which shows some of the features on the south-western part of the image (I stopped read_bin.m at line 57 after fclose(f) to avoid searching for the strong target for azimuth compression).
I'll try to finish the azimuth compression as mentioned to LucaMancaITA by using the theoretical parabolic shape of the phase along azimuth to see if the image can be completely focused in both directions.
Best, JM
I tried
function [sol,ssol,mychirp]=read_bin(filename, count, l)
% file name, number of complex data to read, line length
% x=read_bin('resultSW10_T1286748881_NQ11924.bin',2000,11924*2);
% imagesc(abs(x(1:10:end,:)));
pkg load signal
f=fopen(filename,'rb');
if (f<0) v=0;
else
fref=37.53472224
% UPDATE THESE VALUES: here taken from S1A_IW_RAW__0SDV_20210410T055133_20210410T055205_037384_0467D7_5D9B
IW=str2num(filename(9:10)) % swath number
if (IW==10)
fs=3/7*4*fref % case 10 in Range Decimation, p.35 of Packet Protocol Data Unit
TXPRRcode=1605; % upchirp
TXPSFcode=-12335; % start
TXPLcode=1967;
PRI=21859
else if (IW==11)
fs=4/11*4*fref % case 11 in Range Decimation, p.35 of Packet Protocol Data Unit
TXPRRcode=1160; % upchirp
TXPSFcode=-10546; % start
TXPLcode=2327;
PRI=25857
else if (IW==12)
fs=5/16*4*fref % case 12 in Range Decimation, p.35 of Packet Protocol Data Unit
TXPRRcode=1193; % upchirp
TXPSFcode=-9341; % start
TXPLcode=2004;
PRI=22265
else printf("Wrong swath number");
end
end
end
% range compression
TXPRR=TXPRRcode*fref^2/2^21 % MHz/us
TXPSF=TXPRR/4/fref+TXPSFcode/2^14*fref % MHz
TXPL=TXPLcode/fref % us
N=TXPL*fs
tim=linspace(-TXPL/2,TXPL/2,N);
phi1=TXPSF+TXPRR*TXPL/2
phi2=TXPRR/2
nomchip=exp(j*2*pi*(phi1*tim+phi2*tim.^2));
sol=zeros(l,count); % size(t)(2),size(t)(1)); % avoid dynamic memory allocation
for k=1:count % size(t)(1) % 6592 = sweep along azimuth for range compression
t=fread(f,2*l,'float'); % dynamically read to avoid filling memory with a huge t
t=t(1:2:end)+t(2:2:end)*i; % see GNU Radio's read_complex_binary
k
tmp=xcorr(nomchip,t);
tmp=tmp(floor(length(nomchip)/2)+1:length(t)+floor(length(nomchip)/2));
sol(:,k)=tmp;
end
% range compression completed
clear t
fclose(f)
%
%https://celestrak.com/satcat/tle.php?CATNR=39634
%SENTINEL-1A
%1 39634U 14016A 21079.09896275 -.00000050 00000-0 -79699-6 0 9995
%2 39634 98.1819 88.1319 0001317 80.0185 280.1166 14.59198711370760
% https://sentinels.copernicus.eu/web/sentinel/technical-guides/sentinel-1-sar/products-algorithms/level-1/single-look-complex/stripmap
% range resolutio = 3.6 m
% NQ=9975*2*3.1 => ans = 61845 m @ 43.8 mean angle
PRF=fref*1e6/PRI
t=[0:800]/PRF; % time TODO ADJUST
t=t-t(floor(length(t)/2));
l=300/5405; % lambda wavelength
REarth=6371 % km
Hleo=693 % km
Hgeo=35786 % km
period=sqrt((Hleo+REarth)^3/(Hgeo+REarth)^3)*24*3600 % s
orbit=2*pi*(Hleo+REarth)*1000 % m
% TLE: 24/14.5919871*3600=5921.1 s !
Vs=orbit/period % m/s
Ve=Vs/sqrt(1+Hleo/REarth) % m/s
for th=45 % [42.05 45.49 43.8]
theta=th*pi/180; % S6 angle (rad)
R0=Hleo/cos(theta)*1000; % m
Rpp=Ve^2/R0 % m.s^-2
-4*pi/l*Rpp*(t(2)-t(1))
mychirp=-4*pi/l*Rpp/2*t.^2;mychirp=exp(j*mychirp); % phase evolution of azimuth chirp
figure
plot(t,mychirp)
% azimuth compression
ssol=sol'; % avoid dynamic memory allocation
for k=1:l % size(sol)(1) % 19950 = sweep range for azimuth compression
k
tmp=xcorr(mychirp,sol(k,:));
tmp=tmp(floor(length(mychirp)/2)+1:size(sol)(2)+floor(length(mychirp)/2));
ssol(:,k)=tmp;
end
end
endand running with [x,xx,c]=read_bin('resultSW10_T1286748884_NQ11924.bin',1400,11924*2);
I cannot see much difference with the range compressed data imagesc(abs(x(1:10:end,:))',[0 1000]) and the range+azimuth compressed datasets imagesc(abs(xx(:,1:10:end)),[0 1000]) either due to the lack of structure to assess the data quality of because my azimuth chirp is wrong.
i implemented the code as you guided uptil range compression, I was successful in my amazon dataset, I then also implemented it on your Sao Polo data set, as you guided that water resources can help beginners , I got the following results
for this dataset,
resultSW11_T1297499458_NQ11995.bin
which I think is pretty good, as it does show, that some features are present. Thank you so much for your guidance, I`ll now try to move on to azimuth compression, looking forward for your guidance again.
Thank You,
Best Regards
How are you calculating TXPRRcode, TXPLcode and TXPSFcode for different modes?
Thank you,
Best Regards
See Sentinel-1-SAR-Space-Packet-Protocol-Data-Unit.pdf in the doc/ directory, pp.37-39, with the TXPRRcode, TXPLcode and TXPSFcode values extracted from the raw level0 files telemetry. These values are given on the console output of read_file, e.g.
0c1c: 1(1) 65(65) 12(12) c53f: 3(3) Count=24 Len=16326(61..65533)
Time: 1286748881:20918 352ef853(352EF853) WordIndex=b WordVal=1f67 00009b98 0000a6bb BAQ=0c(c) BlockLen=1f(1F)
Decim=9 TXPRR=^84a9=1193 TXPSF=-0x247d=9341 TXPL=000007d4=2004 PRI=000056f9=22265 Polar=7 Typ=0(0) Swath=c NQ=9884
, finished processing 16264 echo
How do you Plot the uncompressed data?
I am trying the following;
f=fopen(resultSW11_T1297499460_NQ11995.bin,'rb');
t=fread(f,2*l,'float');
imagesc((abs(t(1:end,:))),[0 1000])
I just get a solid colored blue result.
More-so, how were you able to plot the SaoPolo_SM_Level 0 image with such good features,
The read_bin file azimuth compression code i am implementing does not provide with good compression.
Thank you,
Best Regards
I have updated the script as I computed the phase but forgot to create a complex azimuth chirp exp(j*phase). The result is not yet good though.
When you imagesc() you must understand that the last argument [0 1000] is the color scale range: leave it empty at first unless you want to zoom on some feature in the color range by for e.g. saturating a strong target.
Hi jmfriedt,
try to use as azimuth reference chirp the following expression:
As first approximation I set the doppler centroid frequency to zero.
With this I was able to obtain this SLC product (which still it's a bit blurry if compared to the one obtained by using the chirp found in the range-compressed image)
Hope this can be helpful,
Luca
Hi Luca,
Can you share your M.file or code please
Thank you,
best Regards
I updated https://github.com/jmfriedt/sentinel1_level0/tree/main/saopaolo_SM with the script for range compression, extracting the azimuth chirp from the strongest target and azimuth compression. Will try now with this expression of the azimuth chirp. Thanks for the contribution.
I don't know what I did wrong previously, I suspect using exp(phi) instead of exp(j*phi) when writing the phase expression. I have updated https://github.com/jmfriedt/sentinel1_level0/tree/main/saopaolo_SM including pulse_compression.m and the resulting azimuth compressed images in the associated README.md file. Unfortunately I cannot commit all1500bot.mat which is 1.5 GB and even compressed larger than 1 GB: it is the subset of the Sao Paulo Stripmap IQ stream reorganized as a matrix according to the Pulse Repetition Interval, whos on the sol variable:
c sol 19950x1500 478800000 double