Download this notebook from github.

what is specific to TOPS SLC products?

Contrarily to GRD products SLC ones preserve the signal described in time per burst.

It means that a portion of SAR image can be seen up to 4 times by the sensor.

The different bursts are overlapping, and the subswaths are also overlapping.

SLC products also contains the phase information, and digital number in .tiff products are complex values.

How to read a TOPS SLC product with xsar ?

the TOPS (IW and EW) SLC products are distributed by ESA in .SAFE directories containing measurement and annotations per subswath and polarizations.

to open a multi dataset

[1]:

import xsar
import geoviews as gv
import holoviews as hv
import geoviews.feature as gf
hv.extension('bokeh')
path = xsar.get_test_file('S1A_IW_SLC__1SDV_20170907T102951_20170907T103021_018268_01EB76_Z020.SAFE')
sub_swath_IDs = xsar.Sentinel1Meta(path).subdatasets
sub_swath_IDs

[1]:
geometry
SENTINEL1_DS:/tmp/S1A_IW_SLC__1SDV_20170907T102951_20170907T103021_018268_01EB76_Z020.SAFE:IW1 POLYGON ((-67.47144 22.48036, -68.29402 22.627...
SENTINEL1_DS:/tmp/S1A_IW_SLC__1SDV_20170907T102951_20170907T103021_018268_01EB76_Z020.SAFE:IW2 POLYGON ((-68.28370 22.56667, -69.14890 22.716...
SENTINEL1_DS:/tmp/S1A_IW_SLC__1SDV_20170907T102951_20170907T103021_018268_01EB76_Z020.SAFE:IW3 POLYGON ((-69.13847 22.65528, -69.90347 22.784... 
[2]:

multids = xsar.Sentinel1Meta(path)
#print(dir(multids.subdatasets))
for subswath in multids.subdatasets.index:
    print(subswath)
    onesubswath = xsar.Sentinel1Dataset(subswath)
    print(onesubswath.get_bursts_polygons()['geometry'])

SENTINEL1_DS:/tmp/S1A_IW_SLC__1SDV_20170907T102951_20170907T103021_018268_01EB76_Z020.SAFE:IW1

/home3/homedir7/perso/sarwing/xsar/xsar/src/xsar/base_dataset.py:410: FutureWarning: You are adding a column named 'geometry' to a GeoDataFrame constructed without an active geometry column. Currently, this automatically sets the active geometry column to 'geometry' but in the future that will no longer happen. Instead, either provide geometry to the GeoDataFrame constructor (GeoDataFrame(... geometry=GeoSeries()) or use `set_geometry('geometry')` to explicitly set the active geometry column.
  blocks['geometry'] = blocks['geometry_image'].apply(self.coords2ll)

burst
0    POLYGON ((-67.51030 22.31627, -68.32788 22.462...
1    POLYGON ((-67.54546 22.15003, -68.36202 22.296...
2    POLYGON ((-67.58064 21.98355, -68.39618 22.130...
3    POLYGON ((-67.61584 21.81681, -68.43038 21.963...
4    POLYGON ((-67.65099 21.65019, -68.46453 21.797...
5    POLYGON ((-67.68614 21.48344, -68.49869 21.630...
6    POLYGON ((-67.71973 21.31678, -68.53187 21.463...
7    POLYGON ((-67.75480 21.15015, -68.56597 21.297...
8    POLYGON ((-67.78979 20.98376, -68.59997 21.131...
9    POLYGON ((-67.82474 20.81746, -68.63397 20.964...
Name: geometry, dtype: geometry
SENTINEL1_DS:/tmp/S1A_IW_SLC__1SDV_20170907T102951_20170907T103021_018268_01EB76_Z020.SAFE:IW2

/home3/homedir7/perso/sarwing/xsar/xsar/src/xsar/base_dataset.py:410: FutureWarning: You are adding a column named 'geometry' to a GeoDataFrame constructed without an active geometry column. Currently, this automatically sets the active geometry column to 'geometry' but in the future that will no longer happen. Instead, either provide geometry to the GeoDataFrame constructor (GeoDataFrame(... geometry=GeoSeries()) or use `set_geometry('geometry')` to explicitly set the active geometry column.
  blocks['geometry'] = blocks['geometry_image'].apply(self.coords2ll)

burst
0    POLYGON ((-68.32084 22.40319, -69.18158 22.552...
1    POLYGON ((-68.35502 22.23689, -69.21468 22.386...
2    POLYGON ((-68.38922 22.07036, -69.24782 22.220...
3    POLYGON ((-68.42341 21.90382, -69.28095 22.053...
4    POLYGON ((-68.45756 21.73740, -69.31405 21.887...
5    POLYGON ((-68.49178 21.57049, -69.34724 21.720...
6    POLYGON ((-68.52465 21.40395, -69.37954 21.554...
7    POLYGON ((-68.55877 21.23740, -69.41263 21.388...
8    POLYGON ((-68.59282 21.07108, -69.44568 21.221...
9    POLYGON ((-68.62691 20.90451, -69.47877 21.055...
Name: geometry, dtype: geometry
SENTINEL1_DS:/tmp/S1A_IW_SLC__1SDV_20170907T102951_20170907T103021_018268_01EB76_Z020.SAFE:IW3
burst
0    POLYGON ((-69.17395 22.49225, -69.93513 22.620...
1    POLYGON ((-69.20712 22.32585, -69.96735 22.454...
2    POLYGON ((-69.24028 22.15942, -69.99957 22.288...
3    POLYGON ((-69.27343 21.99300, -70.03178 22.121...
4    POLYGON ((-69.30655 21.82670, -70.06397 21.955...
5    POLYGON ((-69.33978 21.65979, -70.09628 21.788...
6    POLYGON ((-69.37186 21.49329, -70.12790 21.622...
7    POLYGON ((-69.40496 21.32697, -70.16009 21.456...
8    POLYGON ((-69.43808 21.16053, -70.19231 21.290...
9    POLYGON ((-69.47121 20.99397, -70.22456 21.123...
Name: geometry, dtype: geometry

/home3/homedir7/perso/sarwing/xsar/xsar/src/xsar/base_dataset.py:410: FutureWarning: You are adding a column named 'geometry' to a GeoDataFrame constructed without an active geometry column. Currently, this automatically sets the active geometry column to 'geometry' but in the future that will no longer happen. Instead, either provide geometry to the GeoDataFrame constructor (GeoDataFrame(... geometry=GeoSeries()) or use `set_geometry('geometry')` to explicitly set the active geometry column.
  blocks['geometry'] = blocks['geometry_image'].apply(self.coords2ll)

[3]:

#gv.tile_sources.Wikipedia*gv.Polygons(multids.bursts(only_valid_location=True)['geometry']).opts(width=800,height=450,alpha=0.5)
import pandas as pd
tmp = pd.concat([xsar.Sentinel1Dataset(onesubswath).get_bursts_polygons(only_valid_location=True) for onesubswath in multids.subdatasets.index])
gv.tile_sources.Wikipedia*gv.Polygons(tmp['geometry']).opts(width=800,height=450,alpha=0.5)

/home3/homedir7/perso/sarwing/xsar/xsar/src/xsar/base_dataset.py:410: FutureWarning: You are adding a column named 'geometry' to a GeoDataFrame constructed without an active geometry column. Currently, this automatically sets the active geometry column to 'geometry' but in the future that will no longer happen. Instead, either provide geometry to the GeoDataFrame constructor (GeoDataFrame(... geometry=GeoSeries()) or use `set_geometry('geometry')` to explicitly set the active geometry column.
  blocks['geometry'] = blocks['geometry_image'].apply(self.coords2ll)
/home3/homedir7/perso/sarwing/xsar/xsar/src/xsar/base_dataset.py:410: FutureWarning: You are adding a column named 'geometry' to a GeoDataFrame constructed without an active geometry column. Currently, this automatically sets the active geometry column to 'geometry' but in the future that will no longer happen. Instead, either provide geometry to the GeoDataFrame constructor (GeoDataFrame(... geometry=GeoSeries()) or use `set_geometry('geometry')` to explicitly set the active geometry column.
  blocks['geometry'] = blocks['geometry_image'].apply(self.coords2ll)
/home3/homedir7/perso/sarwing/xsar/xsar/src/xsar/base_dataset.py:410: FutureWarning: You are adding a column named 'geometry' to a GeoDataFrame constructed without an active geometry column. Currently, this automatically sets the active geometry column to 'geometry' but in the future that will no longer happen. Instead, either provide geometry to the GeoDataFrame constructor (GeoDataFrame(... geometry=GeoSeries()) or use `set_geometry('geometry')` to explicitly set the active geometry column.
  blocks['geometry'] = blocks['geometry_image'].apply(self.coords2ll)
/tmp/ipykernel_25243/3208064731.py:4: GeoviewsDeprecationWarning: 'Wikipedia' is deprecated and will be removed in version 1.11, use 'OSM' instead.
  gv.tile_sources.Wikipedia*gv.Polygons(tmp['geometry']).opts(width=800,height=450,alpha=0.5)

[3]:

[4]:

tmp2 = pd.concat([xsar.Sentinel1Dataset(onesubswath).get_bursts_polygons(only_valid_location=False) for onesubswath in multids.subdatasets.index])
gv.tile_sources.Wikipedia*gv.Polygons(tmp2['geometry']).opts(width=800,height=450,alpha=0.2)*gv.Polygons(tmp['geometry']).opts(width=800,height=450,alpha=0.2)

/home3/homedir7/perso/sarwing/xsar/xsar/src/xsar/base_dataset.py:410: FutureWarning: You are adding a column named 'geometry' to a GeoDataFrame constructed without an active geometry column. Currently, this automatically sets the active geometry column to 'geometry' but in the future that will no longer happen. Instead, either provide geometry to the GeoDataFrame constructor (GeoDataFrame(... geometry=GeoSeries()) or use `set_geometry('geometry')` to explicitly set the active geometry column.
  blocks['geometry'] = blocks['geometry_image'].apply(self.coords2ll)
/home3/homedir7/perso/sarwing/xsar/xsar/src/xsar/base_dataset.py:410: FutureWarning: You are adding a column named 'geometry' to a GeoDataFrame constructed without an active geometry column. Currently, this automatically sets the active geometry column to 'geometry' but in the future that will no longer happen. Instead, either provide geometry to the GeoDataFrame constructor (GeoDataFrame(... geometry=GeoSeries()) or use `set_geometry('geometry')` to explicitly set the active geometry column.
  blocks['geometry'] = blocks['geometry_image'].apply(self.coords2ll)
/home3/homedir7/perso/sarwing/xsar/xsar/src/xsar/base_dataset.py:410: FutureWarning: You are adding a column named 'geometry' to a GeoDataFrame constructed without an active geometry column. Currently, this automatically sets the active geometry column to 'geometry' but in the future that will no longer happen. Instead, either provide geometry to the GeoDataFrame constructor (GeoDataFrame(... geometry=GeoSeries()) or use `set_geometry('geometry')` to explicitly set the active geometry column.
  blocks['geometry'] = blocks['geometry_image'].apply(self.coords2ll)
/tmp/ipykernel_25243/2354858605.py:2: GeoviewsDeprecationWarning: 'Wikipedia' is deprecated and will be removed in version 1.11, use 'OSM' instead.
  gv.tile_sources.Wikipedia*gv.Polygons(tmp2['geometry']).opts(width=800,height=450,alpha=0.2)*gv.Polygons(tmp['geometry']).opts(width=800,height=450,alpha=0.2)

[4]:

to open a specific subswath

to open a subswath one can use xsar.Sentinel1Meta class to get an overview of the meta-data.

[5]:

s1meta = xsar.sentinel1_meta.Sentinel1Meta("SENTINEL1_DS:%s:IW3" % path)
s1meta

[5]:

<Sentinel1Meta single object>

to manipulate the data, the user have to open Sentinel1Dataset instance.

[6]:

s1ds = xsar.sentinel1_dataset.Sentinel1Dataset("SENTINEL1_DS:%s:IW3" % path)
s1ds

[6]:

<Sentinel1Dataset full coverage object>

content of the Sentinel1Dataset python object

  • a datatree (xarray)

  • a dataset (xarray)

  • a s1meta object (class instance)

open the datatree

[7]:

s1ds.datatree

[7]:

<xarray.DatasetView>
Dimensions:  ()
Data variables:
    *empty*
Attributes: (12/15)
    name:              SENTINEL1_DS:/tmp/S1A_IW_SLC__1SDV_20170907T102951_201...
    short_name:        SENTINEL1_DS:S1A_IW_SLC__1SDV_20170907T102951_20170907...
    product:           SLC
    safe:              S1A_IW_SLC__1SDV_20170907T102951_20170907T103021_01826...
    swath:             IW
    multidataset:      False
    ...                ...
    start_date:        2017-09-07 10:29:53.785096
    stop_date:         2017-09-07 10:30:21.707772
    footprint:         POLYGON ((-69.13847173718773 22.65528072350044, -69.90...
    coverage:          187km * 79km (line * sample )
    orbit_pass:        Descending
    platform_heading:  -167.7216262988962

open dataset

[8]:

s1ds.dataset

[8]:

<xarray.Dataset>
Dimensions:         (pol: 2, line: 15070, sample: 23768)
Coordinates:
  * pol             (pol) object 'VV' 'VH'
  * line            (line) int64 0 1 2 3 4 5 ... 15065 15066 15067 15068 15069
  * sample          (sample) int64 0 1 2 3 4 5 ... 23763 23764 23765 23766 23767
Data variables:
    digital_number  (pol, line, sample) complex64 dask.array<chunksize=(1, 5000, 5000), meta=np.ndarray>
    time            (line) datetime64[ns] 2017-09-07T10:29:53.785278 ... 2017...
    sampleSpacing   float64 2.33
    lineSpacing     float64 13.96
Attributes: (12/15)
    name:              SENTINEL1_DS:/tmp/S1A_IW_SLC__1SDV_20170907T102951_201...
    short_name:        SENTINEL1_DS:S1A_IW_SLC__1SDV_20170907T102951_20170907...
    product:           SLC
    safe:              S1A_IW_SLC__1SDV_20170907T102951_20170907T103021_01826...
    swath:             IW
    multidataset:      False
    ...                ...
    start_date:        2017-09-07 10:29:53.785096
    stop_date:         2017-09-07 10:30:21.707772
    footprint:         POLYGON ((-69.13847173718773 22.65528072350044, -69.90...
    coverage:          187km * 79km (line * sample )
    orbit_pass:        Descending
    platform_heading:  -167.7216262988962

get the sar_meta

[9]:

s1ds.sar_meta

[9]:

<BlockingActorProxy: Sentinel1Meta>

add high resolution interpolated variables

for instance longitudes, latitudes, incidence angle,…

Note that this step is automatically done by default for GRD products but not for SLC.

[10]:

s1ds.add_high_resolution_variables()
s1ds.dataset

/home3/homedir7/perso/sarwing/xsar/xsar/src/xsar/base_dataset.py:410: FutureWarning: You are adding a column named 'geometry' to a GeoDataFrame constructed without an active geometry column. Currently, this automatically sets the active geometry column to 'geometry' but in the future that will no longer happen. Instead, either provide geometry to the GeoDataFrame constructor (GeoDataFrame(... geometry=GeoSeries()) or use `set_geometry('geometry')` to explicitly set the active geometry column.
  blocks['geometry'] = blocks['geometry_image'].apply(self.coords2ll)

[10]:

<xarray.Dataset>
Dimensions:               (pol: 2, line: 15070, sample: 23768)
Coordinates:
  * pol                   (pol) object 'VV' 'VH'
  * line                  (line) int64 0 1 2 3 4 ... 15066 15067 15068 15069
  * sample                (sample) int64 0 1 2 3 4 ... 23764 23765 23766 23767
    spatial_ref           int64 0
Data variables: (12/20)
    digital_number        (pol, line, sample) complex64 dask.array<chunksize=(1, 5000, 5000), meta=np.ndarray>
    time                  (line) datetime64[ns] 2017-09-07T10:29:53.785278 .....
    sampleSpacing         float64 2.33
    lineSpacing           float64 13.96
    gamma0_lut            (pol, line, sample) float64 dask.array<chunksize=(1, 4096, 4096), meta=np.ndarray>
    noise_lut_azi         (pol, line, sample) float32 dask.array<chunksize=(1, 4096, 4096), meta=np.ndarray>
    ...                    ...
    incidence             (line, sample) float64 dask.array<chunksize=(4096, 4096), meta=np.ndarray>
    elevation             (line, sample) float64 dask.array<chunksize=(4096, 4096), meta=np.ndarray>
    longitude             (line, sample) float64 dask.array<chunksize=(4096, 4096), meta=np.ndarray>
    latitude              (line, sample) float64 dask.array<chunksize=(4096, 4096), meta=np.ndarray>
    velocity              (line) float64 7.597e+03 7.597e+03 ... 7.597e+03
    range_ground_spacing  (sample) float64 dask.array<chunksize=(4096,), meta=np.ndarray>
Attributes: (12/15)
    name:              SENTINEL1_DS:/tmp/S1A_IW_SLC__1SDV_20170907T102951_201...
    short_name:        SENTINEL1_DS:S1A_IW_SLC__1SDV_20170907T102951_20170907...
    product:           SLC
    safe:              S1A_IW_SLC__1SDV_20170907T102951_20170907T103021_01826...
    swath:             IW
    multidataset:      False
    ...                ...
    start_date:        2017-09-07 10:29:53.785096
    stop_date:         2017-09-07 10:30:21.707772
    footprint:         POLYGON ((-69.13847173718773 22.65528072350044, -69.90...
    coverage:          187km * 79km (line * sample )
    orbit_pass:        Descending
    platform_heading:  -167.7216262988962

alternatively compute the value of a geolocation field (longitude, latitude, incidence,..) at given image coordinates

this method has been added to avoid adding hig resolution grids for some variable while we only need specific points

[11]:

s1ds.get_ll_from_SLC_geoloc(line=5,sample=3000,varname='incidence')

[11]:

array(42.28017613)

perform calibration and denoizing on sigma0, gamma0 and beta0

[12]:

s1ds.apply_calibration_and_denoising()
s1ds.dataset

[12]:

<xarray.Dataset>
Dimensions:               (pol: 2, line: 15070, sample: 23768)
Coordinates:
  * pol                   (pol) object 'VV' 'VH'
  * line                  (line) int64 0 1 2 3 4 ... 15066 15067 15068 15069
  * sample                (sample) int64 0 1 2 3 4 ... 23764 23765 23766 23767
    spatial_ref           int64 0
Data variables: (12/26)
    digital_number        (pol, line, sample) complex64 dask.array<chunksize=(1, 5000, 5000), meta=np.ndarray>
    time                  (line) datetime64[ns] 2017-09-07T10:29:53.785278 .....
    sampleSpacing         float64 2.33
    lineSpacing           float64 13.96
    gamma0_lut            (pol, line, sample) float64 dask.array<chunksize=(1, 4096, 4096), meta=np.ndarray>
    noise_lut_azi         (pol, line, sample) float32 dask.array<chunksize=(1, 4096, 4096), meta=np.ndarray>
    ...                    ...
    sigma0_raw            (pol, line, sample) float64 dask.array<chunksize=(1, 4096, 4096), meta=np.ndarray>
    nesz                  (pol, line, sample) float64 dask.array<chunksize=(1, 4096, 4096), meta=np.ndarray>
    gamma0_raw            (pol, line, sample) float64 dask.array<chunksize=(1, 4096, 4096), meta=np.ndarray>
    negz                  (pol, line, sample) float64 dask.array<chunksize=(1, 4096, 4096), meta=np.ndarray>
    sigma0                (pol, line, sample) float64 dask.array<chunksize=(1, 4096, 4096), meta=np.ndarray>
    gamma0                (pol, line, sample) float64 dask.array<chunksize=(1, 4096, 4096), meta=np.ndarray>
Attributes: (12/15)
    name:              SENTINEL1_DS:/tmp/S1A_IW_SLC__1SDV_20170907T102951_201...
    short_name:        SENTINEL1_DS:S1A_IW_SLC__1SDV_20170907T102951_20170907...
    product:           SLC
    safe:              S1A_IW_SLC__1SDV_20170907T102951_20170907T103021_01826...
    swath:             IW
    multidataset:      False
    ...                ...
    start_date:        2017-09-07 10:29:53.785096
    stop_date:         2017-09-07 10:30:21.707772
    footprint:         POLYGON ((-69.13847173718773 22.65528072350044, -69.90...
    coverage:          187km * 79km (line * sample )
    orbit_pass:        Descending
    platform_heading:  -167.7216262988962
[13]:

s1ds.dataset['ground_heading']

[13]:

<xarray.DataArray 'ground_heading' (line: 15070, sample: 23768)>
dask.array<ground_heading, shape=(15070, 23768), dtype=float32, chunksize=(4096, 4096), chunktype=numpy.ndarray>
Coordinates:
  * line         (line) int64 0 1 2 3 4 5 ... 15065 15066 15067 15068 15069
  * sample       (sample) int64 0 1 2 3 4 5 ... 23763 23764 23765 23766 23767
    spatial_ref  int64 0
Attributes:
    comment:    at ground level, computed from lon/lat in azimuth direction
    long_name:  Platform heading (azimuth from North)
    units:      Degrees

get azimuth time variable

the variable azimuth_time is given in annotations .xml files, it describes the date of acquisition of each pixel in the subswath. it is also called SAR “long time”, in opposition to the “short time” given by the slant_range_time variable

[14]:

aziHr = s1ds.dataset['time']
hv.Curve(aziHr.values)

[14]:

azimuth time variable estimated at the middle (in range) of the dataset selected is showing the bursts overlapping.

This variable is used to rasterize the 7 variables (longitude, latitude,…) described in the geolocationGrid at low resolution.

get the ground range spacing

One specificity of the SLC products is that the ground range spacing is not equal to the slant range spacing (it is the case in GRD products).

[15]:

#print(s1ds.sar_meta.image['ground_pixel_spacing'])
#print(s1ds.sar_meta.image['slant_pixel_spacing'])
print(s1ds.dataset['sampleSpacing'])
print(s1ds.dataset['lineSpacing'])

<xarray.DataArray 'sampleSpacing' ()>
array(2.329562)
Coordinates:
    spatial_ref  int64 0
Attributes:
    units:        m
    referential:  slant
<xarray.DataArray 'lineSpacing' ()>
array(13.95527)
Coordinates:
    spatial_ref  int64 0
Attributes:
    units:    m

The ground range spacing depends of the incidence angle.

\[grdRangeSpacing = \frac{slantRangeSpacing}{sinus(\theta)}\]

It is possible for the users to get the ground range spacing vector along the range axis.

[16]:

rgs = s1ds.dataset['range_ground_spacing']
rgs

[16]:

<xarray.DataArray 'range_ground_spacing' (sample: 23768)>
dask.array<divide, shape=(23768,), dtype=float64, chunksize=(4096,), chunktype=numpy.ndarray>
Coordinates:
  * sample       (sample) int64 0 1 2 3 4 5 ... 23763 23764 23765 23766 23767
    spatial_ref  int64 0
Attributes:
    history:  
[17]:

hv.Curve(rgs).opts(width=400,show_grid=True)

[17]:

get complex digital number

[18]:

s1ds.datatree['measurement']['digital_number']

[18]:

<xarray.DataArray 'digital_number' (pol: 2, line: 15070, sample: 23768)>
dask.array<getitem, shape=(2, 15070, 23768), dtype=complex64, chunksize=(1, 5000, 5000), chunktype=numpy.ndarray>
Coordinates:
  * pol          (pol) object 'VV' 'VH'
  * line         (line) int64 0 1 2 3 4 5 ... 15065 15066 15067 15068 15069
  * sample       (sample) int64 0 1 2 3 4 5 ... 23763 23764 23765 23766 23767
    spatial_ref  int64 0
Attributes:
    comment:  denoised digital number, read at full resolution
    history:  digital_number: measurement/s1a-iw3-slc-v*-20170907t102953-2017...

equivalent to

[19]:

s1ds.dataset['digital_number']

[19]:

<xarray.DataArray 'digital_number' (pol: 2, line: 15070, sample: 23768)>
dask.array<getitem, shape=(2, 15070, 23768), dtype=complex64, chunksize=(1, 5000, 5000), chunktype=numpy.ndarray>
Coordinates:
  * pol          (pol) object 'VV' 'VH'
  * line         (line) int64 0 1 2 3 4 5 ... 15065 15066 15067 15068 15069
  * sample       (sample) int64 0 1 2 3 4 5 ... 23763 23764 23765 23766 23767
    spatial_ref  int64 0
Attributes:
    comment:  denoised digital number, read at full resolution
    history:  digital_number: measurement/s1a-iw3-slc-v*-20170907t102953-2017...

additional informations