Download this notebook from github.


A notebook to illustrate how to compute and correct Impulse Response in IW SLC product

compute IW1 impulse response

[1]:
import warnings
warnings.filterwarnings("ignore")
from xsarslc.processing.impulseResponse import generate_IWS_AUX_file_ImpulseReponse
# SENTINEL1_DS:/home/datawork-cersat-public/cache/project/mpc-sentinel1/data/esa/sentinel-1a/L1/WV/S1A\_WV\_SLC\_\_1S/2018/050/S1A\_WV\_SLC\_\_1SSV\_20180219T221522\_20180219T222851\_020681\_0236CE\_8F55.SAFE:WV_051
subswathes = {'/home/datawork-cersat-public/cache/project/mpc-sentinel1/data/esa/sentinel-1a/L1/IW/S1A_IW_SLC__1S/2022/271/S1A_IW_SLC__1SDV_20220928T095555_20220928T095622_045203_056722_FB4F.SAFE/':[1],
             }
subswath_number = 1

IRs_IW1_VV = generate_IWS_AUX_file_ImpulseReponse(subswathes, subswath_number , polarization='VV')

IRs_IW1_VV
IRs_IW1_VH = generate_IWS_AUX_file_ImpulseReponse(subswathes, subswath_number , polarization='VH')
IRs_IW1_VH
[1]:
<xarray.Dataset> Size: 37kB
Dimensions:         (k_srg: 2048, k_az: 256)
Coordinates:
    pol             <U2 8B 'VH'
  * k_srg           (k_srg) float64 16kB -1.349 -1.347 -1.346 ... 1.346 1.347
  * k_az            (k_az) float64 2kB -0.2236 -0.2219 -0.2201 ... 0.2201 0.2219
Data variables:
    range_IR        (k_srg) float64 16kB 0.0006532 0.000653 ... 0.0006555
    mean_incidence  float64 8B 33.5
    azimuth_IR      (k_az) float64 2kB 0.002633 0.002634 ... 0.002637 0.002631
Attributes: (12/21)
    long_name:              Impulse response in range direction
    name:                   SENTINEL1_DS:/home/datawork-cersat-public/cache/p...
    short_name:             SENTINEL1_DS:S1A_IW_SLC__1SDV_20220928T095555_202...
    product:                SLC
    safe:                   S1A_IW_SLC__1SDV_20220928T095555_20220928T095622_...
    swath:                  IW
    ...                     ...
    platform_heading:       -168.0316299165244
    comment:                denoised digital number, read at full resolution
    history:                digital_number: measurement/s1a-iw1-slc-v*-202209...
    radar_frequency:        5405000454.33435
    azimuth_time_interval:  0.002055556299999998
    subswath:               IW1

save IR file to netCDF

[2]:
import os
out_file = "/tmp/IRs/S1_IR_IW1_VV.nc"
if not os.path.exists("/tmp/IRs/"):
    os.makedirs("/tmp/IRs/")
if 'footprint' in IRs_IW1_VV.attrs:
    IRs_IW1_VV.attrs.update({'footprint': str(IRs_IW1_VV.footprint)})
if 'multidataset' in IRs_IW1_VV.attrs:
    IRs_IW1_VV.attrs.update({'multidataset': str(IRs_IW1_VV.multidataset)})
IRs_IW1_VV.to_netcdf(out_file)
print(out_file)
/tmp/IRs/S1_IR_IW1_VV.nc
[3]:
out_file = "/tmp/IRs/S1_IR_IW1_VH.nc"
if 'footprint' in IRs_IW1_VH.attrs:
    IRs_IW1_VH.attrs.update({'footprint': str(IRs_IW1_VH.footprint)})
if 'multidataset' in IRs_IW1_VH.attrs:
    IRs_IW1_VH.attrs.update({'multidataset': str(IRs_IW1_VH.multidataset)})
IRs_IW1_VH.to_netcdf(out_file)
print(out_file)
/tmp/IRs/S1_IR_IW1_VH.nc

concat IR

[4]:
from xsarslc.processing.impulseResponse import merge_IR_polarization
merge_IR_polarization(IR_directory='/tmp/IRs/')
Merging Impulse Responses done
[4]:
'/tmp/IRs/'

display IW1 range impulse response for VV and VH

[5]:
from matplotlib import pyplot as plt
IRs_IW1_VV['range_IR'].plot(label='VV')
#IRs_IW1_VH['range_IR'].plot(label='VH',alpha=0.7)
plt.legend()
plt.grid(True)
plt.title('subswath = %s'%(IRs_IW1_VV.attrs['subswath']))
plt.show()
../_images/examples_example_IW_compute_and_correct_from_impulse_response_9_0.png

display IW1 azimuth impulse response for VV and VH

[6]:
from matplotlib import pyplot as plt
IRs_IW1_VV['azimuth_IR'].plot(label='VV')
#IRs_IW1_VH['azimuth_IR'].plot(label='VH',alpha=0.7)
plt.legend()
plt.grid(True)
plt.title('subswath = %s'%(IRs_IW1_VV.attrs['subswath']))
plt.show()
../_images/examples_example_IW_compute_and_correct_from_impulse_response_11_0.png

apply IR correction when computing a cross spectra

[7]:
import xsarslc
import xsar
import os
one_tiff = '/home/datawork-cersat-public/cache/project/mpc-sentinel1/data/esa/sentinel-1a/L1/IW/S1A_IW_SLC__1S/2022/127/S1A_IW_SLC__1SDV_20220507T162411_20220507T162439_043107_0525DE_734D.SAFE/measurement/s1a-iw1-slc-vv-20220507t162411-20220507t162439-043107-0525de-004.tiff'
fullpathsafeSLC = os.path.dirname(os.path.dirname(one_tiff))
imagette_number = os.path.basename(one_tiff).split('-')[1].replace('iw','')
str_gdal = 'SENTINEL1_DS:%s:IW%s'%(fullpathsafeSLC,imagette_number)
print('str_gdal',str_gdal)
xsarobj = xsar.Sentinel1Dataset(str_gdal) #,resolution='30m'
xsarobj
str_gdal SENTINEL1_DS:/home/datawork-cersat-public/cache/project/mpc-sentinel1/data/esa/sentinel-1a/L1/IW/S1A_IW_SLC__1S/2022/127/S1A_IW_SLC__1SDV_20220507T162411_20220507T162439_043107_0525DE_734D.SAFE:IW1
[7]:
<Sentinel1Dataset full coverage object>
[8]:
dt = xsarobj.datatree
# dt.load()
dt.load()
[8]:
<xarray.DatasetView> Size: 0B
Dimensions:  ()
Data variables:
    *empty*
Attributes: (12/15)
    name:              SENTINEL1_DS:/home/datawork-cersat-public/cache/projec...
    short_name:        SENTINEL1_DS:S1A_IW_SLC__1SDV_20220507T162411_20220507...
    product:           SLC
    safe:              S1A_IW_SLC__1SDV_20220507T162411_20220507T162439_04310...
    swath:             IW
    multidataset:      False
    ...                ...
    start_date:        2022-05-07 16:24:11.657511
    stop_date:         2022-05-07 16:24:39.551410
    footprint:         POLYGON ((-156.007554522451 22.11305598776757, -156.82...
    coverage:          190km * 86km (line * sample )
    orbit_pass:        Descending
    platform_heading:  -167.7204378813183

compute a first time the cross spectra without impulse response

[9]:
import xsarslc.processing.xspectra
intra_xs_without_IR = xsarslc.processing.xspectra.compute_IW_subswath_intraburst_xspectra(dt,
                                         tile_width={'sample': 5.e3, 'line': 5.e3},
                                         tile_overlap={'sample': -10.e3, 'line': -10.e3},
                                         periodo_width = {'sample': 4000, 'line': 4000},
                                        periodo_overlap = {'sample': 2000, 'line': 2000},
                                        polarization='VV',dev=True)
intra_xs_without_IR
intrabursts: 100%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:05<00:00,  5.67s/it]
[9]:
<xarray.Dataset> Size: 44MB
Dimensions:                   (tile_sample: 6, tile_line: 2, bt_thresh: 5,
                               0tau: 3, freq_line: 57, freq_sample: 446,
                               1tau: 2, 2tau: 1, c_sample: 2, c_line: 2,
                               k_gp: 4, phi_hf: 5, lambda_range_max_macs: 11)
Coordinates:
    pol                       <U2 8B 'VV'
  * bt_thresh                 (bt_thresh) int64 40B 5 50 100 150 200
    k_rg                      (tile_line, tile_sample, freq_sample) float64 43kB ...
    k_az                      (freq_line) float64 456B -0.04409 ... 0.04409
    line                      (tile_line) int64 16B 209 1280
    sample                    (tile_line, tile_sample) int64 96B 1274 ... 19158
    longitude                 (tile_line, tile_sample) float64 96B -156.1 ......
    latitude                  (tile_line, tile_sample) float64 96B 22.1 ... 2...
  * k_gp                      (k_gp) int64 32B 1 2 3 4
  * phi_hf                    (phi_hf) int64 40B 1 2 3 4 5
  * lambda_range_max_macs     (lambda_range_max_macs) float64 88B 25.0 ... 275.0
Dimensions without coordinates: tile_sample, tile_line, 0tau, freq_line,
                                freq_sample, 1tau, 2tau, c_sample, c_line
Data variables: (12/30)
    incidence                 (tile_line, tile_sample) float64 96B 31.21 ... ...
    ground_heading            (tile_line, tile_sample) float64 96B -168.9 ......
    sensing_time              (tile_line, tile_sample) datetime64[ns] 96B 202...
    tile_width_sample         (tile_line, tile_sample) int64 96B 5003 ... 5002
    tile_width_line           (tile_line, tile_sample) int64 96B 4998 ... 4998
    sigma0                    (tile_line, tile_sample) float64 96B 0.1109 ......
    ...                        ...
    burst_corner_longitude    (tile_line, c_sample, c_line) float64 64B -156....
    burst_corner_latitude     (tile_line, c_sample, c_line) float64 64B 22.11...
    cwave_params              (tile_line, tile_sample, k_gp, phi_hf) float64 2kB ...
    macs                      (tile_line, tile_sample, lambda_range_max_macs) complex128 2kB ...
    azimuth_cutoff            (tile_line, tile_sample) float64 96B 265.5 ... ...
    azimuth_cutoff_error      (tile_line, tile_sample) float64 96B 0.06312 .....
Attributes: (12/24)
    name:                          SENTINEL1_DS:/home/datawork-cersat-public/...
    short_name:                    SENTINEL1_DS:S1A_IW_SLC__1SDV_20220507T162...
    product:                       SLC
    safe:                          S1A_IW_SLC__1SDV_20220507T162411_20220507T...
    swath:                         IW
    multidataset:                  False
    ...                            ...
    radar_frequency:               5405000454.33435
    azimuth_time_interval:         0.002055556299999998
    required_tile_width_sample:    5000.0
    required_tile_width_line:      5000.0
    required_tile_overlap_sample:  -10000.0
    required_tile_overlap_line:    -10000.0
[10]:
import xarray as xr
IR_path = '/tmp/IRs/S1_IR_IW1.nc' # this is the output file of the method merge_IR_polarization()
new_ir = xr.open_dataset(IR_path,engine='h5netcdf')
new_ir
[10]:
<xarray.Dataset> Size: 55kB
Dimensions:         (pol: 2, k_srg: 2048, k_az: 256)
Coordinates:
  * pol             (pol) <U2 16B 'VH' 'VV'
  * k_srg           (k_srg) float64 16kB -1.349 -1.347 -1.346 ... 1.346 1.347
  * k_az            (k_az) float64 2kB -0.2236 -0.2219 -0.2201 ... 0.2201 0.2219
Data variables:
    range_IR        (pol, k_srg) float64 33kB ...
    mean_incidence  (pol) float64 16B ...
    azimuth_IR      (pol, k_az) float64 4kB ...
Attributes: (12/21)
    long_name:              Impulse response in range direction
    name:                   SENTINEL1_DS:/home/datawork-cersat-public/cache/p...
    short_name:             SENTINEL1_DS:S1A_IW_SLC__1SDV_20220928T095555_202...
    product:                SLC
    safe:                   S1A_IW_SLC__1SDV_20220928T095555_20220928T095622_...
    swath:                  IW
    ...                     ...
    platform_heading:       -168.0316299165244
    comment:                denoised digital number, read at full resolution
    history:                digital_number: measurement/s1a-iw1-slc-v*-202209...
    radar_frequency:        5405000454.33435
    azimuth_time_interval:  0.002055556299999998
    subswath:               IW1
[11]:
xsarobj = xsar.Sentinel1Dataset(str_gdal) #,resolution='30m'
dt = xsarobj.datatree
dt.load()
dt['measurement']
[11]:
<xarray.DatasetView> Size: 5GB
Dimensions:         (pol: 2, line: 14900, sample: 20566)
Coordinates:
  * line            (line) int64 119kB 0 1 2 3 4 ... 14896 14897 14898 14899
  * sample          (sample) int64 165kB 0 1 2 3 4 ... 20562 20563 20564 20565
  * pol             (pol) object 16B 'VV' 'VH'
Data variables:
    digital_number  (pol, line, sample) complex64 5GB 0j 0j 0j 0j ... 0j 0j 0j
    time            (line) datetime64[ns] 119kB 2022-05-07T16:24:11.657340 .....
    sampleSpacing   float64 8B 2.33
    lineSpacing     float64 8B 14.0
Attributes: (12/15)
    name:              SENTINEL1_DS:/home/datawork-cersat-public/cache/projec...
    short_name:        SENTINEL1_DS:S1A_IW_SLC__1SDV_20220507T162411_20220507...
    product:           SLC
    safe:              S1A_IW_SLC__1SDV_20220507T162411_20220507T162439_04310...
    swath:             IW
    multidataset:      False
    ...                ...
    start_date:        2022-05-07 16:24:11.657511
    stop_date:         2022-05-07 16:24:39.551410
    footprint:         POLYGON ((-156.007554522451 22.11305598776757, -156.82...
    coverage:          190km * 86km (line * sample )
    orbit_pass:        Descending
    platform_heading:  -167.7204378813183
[12]:
# print(dt)
polarization = 'VV'
for g in dt.groups:
    # print(g)
    if 'pol' in dt[g].coords:
        print('pol present',g)
        dt[g].sel(pol=polarization)
    else:
        print('pol absent',g)

# 'pol' in dt['intraburst'].coords
pol absent /
pol present /measurement
pol absent /geolocation_annotation
pol absent /bursts
pol absent /FMrate
pol absent /doppler_estimate
pol absent /orbit
pol absent /image
pol present /calibration
pol present /noise_range
pol present /noise_azimuth
pol absent /antenna_pattern
pol absent /swath_merging
pol absent /recalibration

compute a second time the cross spectra with Impulse Response

[13]:
 # compute the same cross spectra but using Impulse Response
intra_xs_with_IR = xsarslc.processing.xspectra.compute_IW_subswath_intraburst_xspectra(dt,
                                         tile_width={'sample': 5.e3, 'line': 5.e3},
                                         tile_overlap={'sample': -10.e3, 'line': -10.e3},
                                         periodo_width = {'sample': 4000, 'line': 4000},
                                        periodo_overlap = {'sample': 2000, 'line': 2000},
                                        polarization='VV',
                                        dev=True,IR_path=IR_path)
intra_xs_with_IR
intrabursts: 100%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:06<00:00,  6.50s/it]
[13]:
<xarray.Dataset> Size: 44MB
Dimensions:                   (tile_sample: 6, tile_line: 2, bt_thresh: 5,
                               0tau: 3, freq_line: 57, freq_sample: 446,
                               1tau: 2, 2tau: 1, c_sample: 2, c_line: 2,
                               k_gp: 4, phi_hf: 5, lambda_range_max_macs: 11)
Coordinates:
    pol                       <U2 8B 'VV'
  * bt_thresh                 (bt_thresh) int64 40B 5 50 100 150 200
    k_rg                      (tile_line, tile_sample, freq_sample) float64 43kB ...
    k_az                      (freq_line) float64 456B -0.04409 ... 0.04409
    line                      (tile_line) int64 16B 209 1280
    sample                    (tile_line, tile_sample) int64 96B 1274 ... 19158
    longitude                 (tile_line, tile_sample) float64 96B -156.1 ......
    latitude                  (tile_line, tile_sample) float64 96B 22.1 ... 2...
  * k_gp                      (k_gp) int64 32B 1 2 3 4
  * phi_hf                    (phi_hf) int64 40B 1 2 3 4 5
  * lambda_range_max_macs     (lambda_range_max_macs) float64 88B 25.0 ... 275.0
Dimensions without coordinates: tile_sample, tile_line, 0tau, freq_line,
                                freq_sample, 1tau, 2tau, c_sample, c_line
Data variables: (12/30)
    incidence                 (tile_line, tile_sample) float64 96B 31.21 ... ...
    ground_heading            (tile_line, tile_sample) float64 96B -168.9 ......
    sensing_time              (tile_line, tile_sample) datetime64[ns] 96B 202...
    tile_width_sample         (tile_line, tile_sample) int64 96B 5003 ... 5002
    tile_width_line           (tile_line, tile_sample) int64 96B 4998 ... 4998
    sigma0                    (tile_line, tile_sample) float64 96B 0.1109 ......
    ...                        ...
    burst_corner_longitude    (tile_line, c_sample, c_line) float64 64B -156....
    burst_corner_latitude     (tile_line, c_sample, c_line) float64 64B 22.11...
    cwave_params              (tile_line, tile_sample, k_gp, phi_hf) float64 2kB ...
    macs                      (tile_line, tile_sample, lambda_range_max_macs) complex128 2kB ...
    azimuth_cutoff            (tile_line, tile_sample) float64 96B 251.7 ... ...
    azimuth_cutoff_error      (tile_line, tile_sample) float64 96B 0.05023 .....
Attributes: (12/24)
    name:                          SENTINEL1_DS:/home/datawork-cersat-public/...
    short_name:                    SENTINEL1_DS:S1A_IW_SLC__1SDV_20220507T162...
    product:                       SLC
    safe:                          S1A_IW_SLC__1SDV_20220507T162411_20220507T...
    swath:                         IW
    multidataset:                  False
    ...                            ...
    radar_frequency:               5405000454.33435
    azimuth_time_interval:         0.002055556299999998
    required_tile_width_sample:    5000.0
    required_tile_width_line:      5000.0
    required_tile_overlap_sample:  -10000.0
    required_tile_overlap_line:    -10000.0
[14]:
tile_line_i = 0
tile_sample_i = 0
onespec_without_IR = intra_xs_without_IR['xspectra_2tau'].isel({'tile_sample':tile_sample_i,'tile_line':tile_line_i}).squeeze()
onespec_without_IR = onespec_without_IR.swap_dims({'freq_line': 'k_az', 'freq_sample': 'k_rg'})
abs(onespec_without_IR.real.mean(dim='k_rg')).plot(label='real without IR')

onespec_with_IR = intra_xs_with_IR['xspectra_2tau'].isel({'tile_sample':tile_sample_i,'tile_line':tile_line_i}).squeeze()
onespec_with_IR = onespec_with_IR.swap_dims({'freq_line': 'k_az', 'freq_sample': 'k_rg'})
abs(onespec_with_IR.real.mean(dim='k_rg')).plot(label='real with IR')
plt.legend()
plt.grid(True)
plt.show()
../_images/examples_example_IW_compute_and_correct_from_impulse_response_22_0.png
[15]:
onespec_without_IR.imag.mean(dim='k_rg').plot(label='imag. without IR')
onespec_with_IR.imag.mean(dim='k_rg').plot(label='imag. with IR')
plt.legend()
plt.grid(True)
plt.show()
../_images/examples_example_IW_compute_and_correct_from_impulse_response_23_0.png