Inputs from Guisepe Grieco (CNR):
You can download it from the SAOCOM web page managed by ASI. Nando, from
the Sapienza university (in copy) should be able to add you to the
licence and give you the credentials. Is that correct Nando?
Concerning the Aquarius GMFs, you can get them from the following link
in the shared ROSS folder:
https://drive.google.com/drive/folders/1_2SdM9NqI5XEJ6RVpAk7vpd8JPGjzK8n?usp=drive_link
[1]
In this folder, you can also find the matlab script to read the files.
Don't hesitate to ask for my support. I guess you will need it.
Cheers,
Giuseppe
SOACOM.py
# -*- coding: utf-8 -*-
"""
Created on Wed Jun 24 15:24:09 2015
@author: imac
"""
import os
import glob
import rasterio as rio
from rasterio.windows import Window
import warnings
import numpy as np
def evalLcCCF(im, dx_az, dx_ra):
''' Returns and estimation of the azimuth wavelength cut-off by means of \n
a least-square fit of the CCF with a Gaussian function \n
Call signature: \n
lc,fxfit,ipsda,kmin,kmax,lmin,lmax=evalLcCCF(im, dl, hp) \n
INPUT: \n
im integer SAR image (electric field) \n
dx_az pixel spacing in azimuth \n
dx_ra pixel spacing in range \n
OUTPUT: \n
lc wavelength cut-off least square method\n
fxfit gaussian fit of the ACF function \n
ipsda ACF \n
kmin minimum wavenumber \n
kmax maximum wavenumber \n
lmin minimum value of the spatial domain \n
lmax masimum value of the spatial domain \n
chi2 chi square \n
pvalue pvalue associated to chi2 \n
Giuseppe Grieco \n
23 Dec 2025
'''
from scipy import optimize
from numpy.fft import ifft, ifft2, fft2, fftshift, ifftshift, fftfreq
#==============================================================================
# Proietto l'immagine nello spazio del Doppler. Definisco
# - numero di look
# - size in azimuth e range
# - costruisco una hanning window per fare un'apodizzazione dello spettro
#==============================================================================
hp_az = int(im.shape[0]/2)
hp_ra = int(im.shape[1]/2)
az_size = 2*hp_az
ran_size = 2*hp_ra
az_win = np.hanning(az_size+2)[1:-1]
ran_win = np.hanning(ran_size+2)[1:-1]
intensity = (abs(im)**2).mean()
window = np.sqrt(1/intensity)
#window = np.sqrt(np.outer(az_win, ran_win)/intensity)
#==============================================================================
# Costruisco l'asse nel dominio della CCF. Definisco
# - pixel spacing in azimuth
#==============================================================================
az_look_size_per = .3
ran_look_size_per = .78
az_look_sep_per = .01
ran_look_sep_per = .02
n_looks_az = int(np.floor((1+az_look_sep_per)/(az_look_size_per+az_look_sep_per)))
n_looks_ran = int(np.floor((1+ran_look_sep_per)/(ran_look_size_per+ran_look_sep_per)))
nsteps = n_looks_az-1
az_look_sep_pt = int(np.round(az_look_sep_per*az_size))
az_look_size_pt = int(az_look_size_per*az_size)
ran_look_size_pt = int(np.round(ran_size*ran_look_size_per))
ran_look_sep_pt = int(np.round(ran_size*ran_look_sep_per))
d_dop = 1/(az_size*dx_az)
x_ccf = fftshift(fftfreq(az_size, d_dop))
image_fft = fftshift(fft2(window*im))
i_az = int((az_size-az_look_size_pt*n_looks_az-az_look_sep_pt*(n_looks_az-1))/2)
i_ran = int(np.ceil((ran_size-ran_look_size_pt)/2))
az_look_iind = np.arange(i_az, az_size, (az_look_size_pt+az_look_sep_pt))[:n_looks_az]
az_look_find = az_look_iind+az_look_size_pt
ran_look_iind = np.arange(i_ran, ran_size, (ran_look_size_pt+ran_look_sep_pt))[:n_looks_ran]
ran_look_find = ran_look_iind+ran_look_size_pt
#==============================================================================
# Inizializzo
# - look nel dominio Doppler
# - Power spectral density
# - Cross Correlation Spectra (CCS)
# - CCS mediati sul numero di step disponibili
# - Cross Correlation Function (CCF)
#==============================================================================
look = np.zeros((az_size, ran_size, n_looks_az, n_looks_ran))
int_look_fft = np.zeros((az_size, ran_size, n_looks_az, n_looks_ran), dtype='complex64')
image_look_fft = np.zeros((az_size, ran_size), dtype='complex64')
int_look = np.zeros((az_size, ran_size, n_looks_az, n_looks_ran), dtype='complex64')
ccs = np.zeros((az_size, ran_size, n_looks_az, n_looks_az), dtype='complex64')
accs = np.zeros((az_size, ran_size, n_looks_az-1), dtype='complex64')
ccf = np.zeros((az_size, nsteps))
#==============================================================================
# Calcolo la densitá spettrale per ogni look
#==============================================================================
for i_az_look in np.arange(n_looks_az):
for i_ran_look in np.arange(n_looks_ran):
image_look_fft = np.zeros((az_size, ran_size), dtype='complex64')
image_look_fft[az_look_iind[i_az_look]:az_look_find[i_az_look], ran_look_iind[i_ran_look]: ran_look_find[i_ran_look]] = image_fft[az_look_iind[i_az_look]:az_look_find[i_az_look], ran_look_iind[i_ran_look]: ran_look_find[i_ran_look]]
look[:, :, i_az_look, i_ran_look] = ifft2(ifftshift(image_look_fft))
dum = fft2(abs(ifft2(image_look_fft))**2)
int_look_fft[:, :, i_az_look, i_ran_look] = dum
int_look[:, :, i_az_look, i_ran_look] = look[:, :, i_az_look, i_ran_look]**2
#==============================================================================
# Calcolo i cross correlation spectra per tutte le combinazioni di look
#==============================================================================
for i_az_look in np.arange(n_looks_az):
for j_az_look in np.arange(i_az_look, n_looks_az):
ccs[:, :, i_az_look, j_az_look] = np.squeeze(int_look_fft[:, :, i_az_look]*np.conj(int_look_fft[:, :, j_az_look]))
#==============================================================================
# Remove DC value
#==============================================================================
ccs[np.tile([0, 1, -1], 3), np.repeat([0, 1, -1], 3), :, :] = 0
int_look_fft[np.tile([0, 1, -1], 3), np.repeat([0, 1, -1], 3), :, :] = 0
for step in np.arange(1, n_looks_az):
for i_az_look in np.arange(0, n_looks_az-step):
accs[:, :, step-1] = accs[:, :, step-1] + 1/(n_looks_az-step)*ccs[:, :, i_az_look, i_az_look+step]
for step in np.arange(nsteps):
dum = fftshift(ifft(ifftshift(accs[:, :, step].mean(axis=1))))
ccf[:, step] = np.real(np.sqrt(dum*np.conj(dum)))
ccf[:, step] = ccf[:, step]/ccf[:, step].max()
p0 = [0., 200/np.sqrt(2)/np.pi, 4., np.mean(np.abs(ccf[:, 0]))]
outFitLS = optimize.leastsq(errfunc, p0, args=(x_ccf, ccf[:, 0]), full_output=1)
lc = np.sqrt(2)*np.pi*outFitLS[0][1]
fxfit = gaussian(outFitLS[0], x_ccf)
#chi2, pvalue = chisquare(np.abs(dum), np.abs(fxfit), ddof=n-4)
chi2 = np.nan
pvalue = np.nan
import pdb; pdb.set_trace()
return lc, fxfit, dum, x_ccf[0], x_ccf[-1], chi2, pvalue, ccf[:, 0], outFitLS[0]
def errfunc(p, x, y):
'''Returns the difference between data and the gaussian function \n
Call signature: \n
d=errfunc(p,x,y) \n
INPUT: \n
p list of the parameters of the Gaussian distribution:\n
p[0]=m mean \n
p[1]=std standard deviation \n
p[2]=amp amplitude factor \n
p[3]=offset offset with respect to the x axis \n
x independent coordinate (space, time, ecc...) \n
y observations \n
OUTPUT: \n
d difference between observations and Gaussian fit: \n
y-{p[3]+p[2]*Gauss(p[0],p[1])} \n
Giuseppe Grieco \n
26 Jun 2015
'''
return y-gaussian(p, x)
def gaussian(B, x):
''' Returns the gaussian function for the given x and the parameters of the distribution\n
Call signature: \n
y=gaussian(B,x) \n
INPUT: \n
B list of the parameters of the distribution \n
m mean of the Gaussian function \n
st standard deviation \n
amp amplification factor \n
offset offest with respect to the x axis. \n
x independent coordinate (space, time, ecc..) \n
OUTPUT: \n
y array of the dependent variable computed as follows \n
y=offset+amp/(st*np.sqrt(2*np.pi))*np.exp(-((x-m)**2/(2*st**2))) \n
Giuseppe Grieco \n
26 Jun 2015
'''
return B[3]+B[2]/(B[1]*np.sqrt(2*np.pi))*np.exp(-((x-B[0])**2/(2*B[1]**2)))
def read_SLC_SM(dir_name, pol='HH', chunk=[0, 1, 0, 1]):
'''
This function reads SAOCOM StripMap SAR images from tiff files.
Note that these images are already radiometrically calibrated.
INPUT:
dir_name name of the directory where the entire SAOCOM file is located
pol polarization ('HH','HV') optional: default='HH'
chunk list containing respectively the first and last indices in the
range direction and in the azimuth direction [iR1,iR2,iA1,iA2].
Default: [0, 1, 0, 1]
OUTPUT:
image subset image from the tiff file
ncols number of columns of the entire image
nrows number of rows of the entire image
'''
#------------------------------------------------------------------------------
# Directory where tiff files are stored in SAOCOM files
#------------------------------------------------------------------------------
meas_dir_name = os.path.join(os.environ['SAO_DATA'], dir_name, 'Data')
tiff_name = glob.glob(os.path.join(meas_dir_name, '*'+pol.lower()))[0]
ras = rio.open(tiff_name)
width = ras.width
height = ras.height
if (chunk[1] > width):
warnings.warn("The final range index is larger than the image width. Exceeding columns will be ignored")
# import pdb; pdb.set_trace()
if (chunk[3] > height):
warnings.warn("The final azimuth index is larger than the image heigth. Exceeding rows will be ignored")
w_size = int(chunk[1]-chunk[0])
h_size = int(chunk[3]-chunk[2])
values = ras.read(1, window=Window(chunk[0], chunk[2], w_size, h_size))
return valuesreadSLC_SOACOM.py
#!/home/giuseppe/.venv/bin/python3
# -*- coding: utf-8 -*-
"""
author: grieco
"""
from SAOCOM import read_SLC_SM, evalLcCCF, errfunc, gaussian
import numpy as np
import os
import matplotlib.pyplot as plt
from scipy import optimize
from numpy.fft import ifft, ifft2, fft2, fftshift, ifftshift, fftfreq
plt.ion()
SAOCOM_DIR = '/home/giuseppe/SAOCOM'
dir_name = os.path.join(SAOCOM_DIR, 'EOL1ASARSAO1A13639337')
pol = 'HH'
dx_az = 5.
dx_ra = 6.
Dx_az = 1e4
Dx_ra = 1e4
dn_az = 2**int(np.round(np.log2(Dx_az/2/dx_az)))
dn_ra = 2**int(np.round(np.log2(Dx_ra/2/dx_ra)))
ic_ra = 3500
ic_az = 9500
chunk = [ic_ra-dn_ra, ic_ra+dn_ra, ic_az-dn_az, ic_az+dn_az]
im = read_SLC_SM(dir_name, pol, chunk)
im_abs = abs(im)
sigma_0 = 10*np.log10(im_abs**2)
# lc, fxfit, dum, x_ccf[0], x_ccf[-1], chi2, pvalue, ccf[:, 0], outFitLS[0] = evalLcCCF(im, dx_az, dx_ra)
#==============================================================================
# Proietto l'immagine nello spazio del Doppler. Definisco
# - numero di look
# - size in azimuth e range
# - costruisco una hanning window per fare un'apodizzazione dello spettro
#==============================================================================
hp_az = int(im.shape[0]/2)
hp_ra = int(im.shape[1]/2)
az_size = 2*hp_az
ran_size = 2*hp_ra
az_win = np.hanning(az_size+2)[1:-1]
ran_win = np.hanning(ran_size+2)[1:-1]
intensity = (im_abs**2).mean()
window = np.sqrt(1/intensity)
# window = np.sqrt(np.outer(az_win, ran_win)/intensity)
#==============================================================================
# Costruisco l'asse nel dominio della CCF. Definisco
# - pixel spacing in azimuth
#==============================================================================
az_look_size_per = .3
ran_look_size_per = .78
az_look_sep_per = .01
ran_look_sep_per = .02
n_looks_az = int(np.floor((1+az_look_sep_per)/(az_look_size_per+az_look_sep_per)))
n_looks_ran = int(np.floor((1+ran_look_sep_per)/(ran_look_size_per+ran_look_sep_per)))
nsteps = n_looks_az-1
az_look_sep_pt = int(np.round(az_look_sep_per*az_size))
az_look_size_pt = int(az_look_size_per*az_size)
ran_look_size_pt = int(np.round(ran_size*ran_look_size_per))
ran_look_sep_pt = int(np.round(ran_size*ran_look_sep_per))
d_dop = 1/(az_size*dx_az)
x_ccf = fftshift(fftfreq(az_size, d_dop))
image_fft = fftshift(fft2(window*im))
i_az = int((az_size-az_look_size_pt*n_looks_az-az_look_sep_pt*(n_looks_az-1))/2)
i_ran = int(np.ceil((ran_size-ran_look_size_pt)/2))
az_look_iind = np.arange(i_az, az_size, (az_look_size_pt+az_look_sep_pt))[:n_looks_az]
az_look_find = az_look_iind+az_look_size_pt
ran_look_iind = np.arange(i_ran, ran_size, (ran_look_size_pt+ran_look_sep_pt))[:n_looks_ran]
ran_look_find = ran_look_iind+ran_look_size_pt
#==============================================================================
# Inizializzo
# - look nel dominio Doppler
# - Power spectral density
# - Cross Correlation Spectra (CCS)
# - CCS mediati sul numero di step disponibili
# - Cross Correlation Function (CCF)
#==============================================================================
look = np.zeros((az_size, ran_size, n_looks_az, n_looks_ran), dtype='complex64')
int_look_fft = np.zeros((az_size, ran_size, n_looks_az, n_looks_ran), dtype='complex64')
image_look_fft = np.zeros((az_size, ran_size), dtype='complex64')
int_look = np.zeros((az_size, ran_size, n_looks_az, n_looks_ran), dtype='complex64')
ccs = np.zeros((az_size, ran_size, n_looks_az, n_looks_az), dtype='complex64')
accs = np.zeros((az_size, ran_size, n_looks_az-1), dtype='complex64')
ccf = np.zeros((az_size, nsteps))
#==============================================================================
# Calcolo la densitá spettrale per ogni look
#==============================================================================
for i_az_look in np.arange(n_looks_az):
for i_ran_look in np.arange(n_looks_ran):
image_look_fft = np.zeros((az_size, ran_size), dtype='complex64')
image_look_fft[az_look_iind[i_az_look]:az_look_find[i_az_look], ran_look_iind[i_ran_look]: ran_look_find[i_ran_look]] = image_fft[az_look_iind[i_az_look]:az_look_find[i_az_look], ran_look_iind[i_ran_look]: ran_look_find[i_ran_look]]
look[:, :, i_az_look, i_ran_look] = ifft2(ifftshift(image_look_fft))
dum = fft2(abs(ifft2(image_look_fft))**2)
int_look_fft[:, :, i_az_look, i_ran_look] = dum
int_look[:, :, i_az_look, i_ran_look] = np.abs(look[:, :, i_az_look, i_ran_look])**2
#==============================================================================
# Calcolo i cross correlation spectra per tutte le combinazioni di look
#==============================================================================
for i_az_look in np.arange(n_looks_az):
for j_az_look in np.arange(i_az_look, n_looks_az):
ccs[:, :, i_az_look, j_az_look] = np.squeeze(int_look_fft[:, :, i_az_look]*np.conj(int_look_fft[:, :, j_az_look]))
#==============================================================================
# Remove DC value
#==============================================================================
ccs[np.tile([0, 1, -1], 3), np.repeat([0, 1, -1], 3), :, :] = 0
int_look_fft[np.tile([0, 1, -1], 3), np.repeat([0, 1, -1], 3), :, :] = 0
for step in np.arange(1, n_looks_az):
for i_az_look in np.arange(0, n_looks_az-step):
accs[:, :, step-1] = accs[:, :, step-1] + 1/(n_looks_az-step)*ccs[:, :, i_az_look, i_az_look+step]
for step in np.arange(nsteps):
dum = fftshift(ifft(ifftshift(accs[:, :, step].mean(axis=1))))
ccf[:, step] = np.real(np.sqrt(dum*np.conj(dum)))
ccf[:, step] = ccf[:, step]/ccf[:, step].max()
p0 = [0., 200/np.sqrt(2)/np.pi, 4., np.mean(np.abs(ccf[:, 0]))]
outFitLS = optimize.leastsq(errfunc, p0, args=(x_ccf, ccf[:, 0]), full_output=1)
lc = np.sqrt(2)*np.pi*outFitLS[0][1]
fxfit = gaussian(outFitLS[0], x_ccf)
#chi2, pvalue = chisquare(np.abs(dum), np.abs(fxfit), ddof=n-4)
chi2 = np.nan
pvalue = np.nan
import pdb; pdb.set_trace()
#==============================================================================
# Plot sigma0
#==============================================================================
vmin = np.percentile(sigma_0, .5)
vmax = np.percentile(sigma_0, 99.5)
plt.figure()
plt.imshow(sigma_0, cmap='gray', vmin=-40., vmax=10.)
#==============================================================================
# Plot abs(image_fft)
#==============================================================================
vmax = np.percentile(np.reshape(abs(image_fft), (image_fft.size,)), 99.5)
lev_dist = vmax/10.
plt.figure()
plt.contour(abs(image_fft), cmap='gray', levels=np.arange(0, vmax+1, lev_dist))
#==============================================================================
# Plot abs averaged cross-spectra for step 0
#==============================================================================
vmax = np.max(abs(accs[:, :, 0]))
lev_dist = vmax/10.
plt.figure()
plt.contour(fftshift(abs(accs[:, :, 0])), cmap='gray', levels=np.arange(0, vmax+1, lev_dist))
#==============================================================================
# Plot abs averaged cross-spectra for step 1
#==============================================================================
vmax = np.max(abs(accs[:, :, 1]))
lev_dist = vmax/10.
plt.figure()
plt.contour(fftshift(abs(accs[:, :, 0])), cmap='gray', levels=np.arange(0, vmax+1, lev_dist))
Inputs from Guisepe Grieco (CNR):
SOACOM.py
readSLC_SOACOM.py