Download this notebook from github.

Advanced XSAR example

[1]:

import xsar
import os
import numpy as np

Sentinel1 example

[2]:

# get test file. You can replace with an path to other SAFE
filename = xsar.get_test_file('S1A_IW_GRDH_1SDV_20170907T103020_20170907T103045_018268_01EB76_Z010.SAFE')

Open a dataset with a xsar.Sentinel1Meta object

A xsar.Sentinel1Meta object handles all attributes and methods that can’t be embdeded in a xarray.Dataset object. It can also replace a filename in xarray.open_dataset

[3]:

sar_meta = xsar.Sentinel1Meta(filename)
sar_meta

[3]:

<Sentinel1Meta single object>

If holoviews extension is loaded, the <Sentinel1Meta objet> have a nice representation. (matplolib is also a valid extension)

[4]:

import holoviews as hv
hv.extension('bokeh')
sar_meta

/home/datawork-cersat-public/cache/project/mpc-sentinel1/workspace/mamba/envs/testgeoviewsnumpy/lib/python3.11/site-packages/dask/dataframe/__init__.py:42: FutureWarning:
Dask dataframe query planning is disabled because dask-expr is not installed.

You can install it with `pip install dask[dataframe]` or `conda install dask`.
This will raise in a future version.

  warnings.warn(msg, FutureWarning)

[4]:

<Sentinel1Meta single object>

sar_meta object is an xsar.Sentinel1Meta object that can be given to xarray.open_dataset or xsar.Sentinel1Dataset , as if it was a filename:

[5]:

sar_ds = xsar.Sentinel1Dataset(sar_meta)
sar_ds

/home/datawork-cersat-public/cache/project/mpc-sentinel1/workspace/mamba/envs/testgeoviewsnumpy/lib/python3.11/site-packages/numpy/_core/numeric.py:452: RuntimeWarning: invalid value encountered in cast
  multiarray.copyto(res, fill_value, casting='unsafe')

[5]:

<Sentinel1Dataset full coverage object>

Open a dataset at lower resolution

resolution keyword can be used to open a dataset at lower resolution.

It might be:

a dict {'line': 20, 'sample': 20}: 20*20 pixels. so if sensor resolution is 10m, the final resolution will be 100m
a string like '200m': Sensor resolution will be automatically used to compute the window size

Then we can instantiate a xsar.Sentinel1Dataset, with the given resolution. Note that the above pixel size has changed.

[6]:

sar_ds = xsar.Sentinel1Dataset(sar_meta, resolution='200m')
sar_ds

/home/datawork-cersat-public/cache/project/mpc-sentinel1/workspace/mamba/envs/testgeoviewsnumpy/lib/python3.11/site-packages/numpy/_core/numeric.py:452: RuntimeWarning: invalid value encountered in cast
  multiarray.copyto(res, fill_value, casting='unsafe')

[6]:

<Sentinel1Dataset full coverage object>

Extract a sub image of 10*10km around a lon/lat point

Convert (lon,lat) to (line, sample)

we can use sar_meta.ll2coords to convert (lon,lat) to (line, sample):

[7]:

# from a shapely object
point_lonlat =  sar_meta.footprint.centroid
point_coords = sar_meta.ll2coords(point_lonlat.x, point_lonlat.y)
point_coords

[7]:

(np.float64(8421.566059451885), np.float64(12589.902956209495))

The result is floating, because it is the position inside the pixel. If real indexes from existing dataset is needed, you’ll have to use sar_ds.ll2coords Result will be the nearest (line, sample) in the dataset

[8]:

point_coords = sar_ds.ll2coords(point_lonlat.x, point_lonlat.y)
point_coords

[8]:

(8429.5, 12589.5)

Extract the sub-image

[9]:

box_size = 10000 # 10km
dist = {'line' : int(np.round(box_size / 2 / sar_meta.pixel_line_m)), 'sample': int(np.round(box_size / 2 / sar_meta.pixel_sample_m))}
dist

[9]:

{'line': 500, 'sample': 500}

The xarray/dask dataset is available as a property : sar_ds.dataset. This attribute can be set to a new values, so the attributes like pixel spacing and coverage are correctly recomputed:

[10]:

# select 10*10 km around point_coords
sar_ds.dataset = sar_ds.dataset.sel(line=slice(point_coords[0] - dist['line'], point_coords[0] + dist['line']), sample=slice(point_coords[1] - dist['sample'], point_coords[1] + dist['sample']))
sar_ds

[10]:

<Sentinel1Dataset sliced object>

[11]:

sar_ds.dataset

[11]:

<xarray.Dataset> Size: 237MB
Dimensions:               (line: 838, sample: 1259, pol: 2)
Coordinates:
  * line                  (line) float64 7kB 9.5 29.5 ... 1.673e+04 1.675e+04
  * sample                (sample) float64 10kB 9.5 29.5 ... 2.515e+04 2.517e+04
  * pol                   (pol) object 16B 'VV' 'VH'
    spatial_ref           int64 8B 0
Data variables: (12/27)
    digital_number        (pol, line, sample) uint16 4MB dask.array<chunksize=(1, 838, 1259), meta=np.ndarray>
    time                  (line) datetime64[ns] 7kB 2017-09-07T10:30:20.95128...
    sampleSpacing         float64 8B 200.0
    lineSpacing           float64 8B 200.0
    gamma0_lut            (pol, line, sample) float64 17MB dask.array<chunksize=(1, 838, 1259), meta=np.ndarray>
    noise_lut_azi         (pol, line, sample) float32 8MB dask.array<chunksize=(1, 838, 1259), meta=np.ndarray>
    ...                    ...
    sigma0_raw            (pol, line, sample) float64 17MB dask.array<chunksize=(1, 838, 1259), meta=np.ndarray>
    nesz                  (pol, line, sample) float64 17MB dask.array<chunksize=(1, 838, 1259), meta=np.ndarray>
    gamma0_raw            (pol, line, sample) float64 17MB dask.array<chunksize=(1, 838, 1259), meta=np.ndarray>
    negz                  (pol, line, sample) float64 17MB dask.array<chunksize=(1, 838, 1259), meta=np.ndarray>
    sigma0                (pol, line, sample) float64 17MB dask.array<chunksize=(1, 838, 1259), meta=np.ndarray>
    gamma0                (pol, line, sample) float64 17MB dask.array<chunksize=(1, 838, 1259), meta=np.ndarray>
Attributes: (12/15)
    name:              SENTINEL1_DS:/home1/scratch/agrouaze/xsardatasync/xsar...
    short_name:        SENTINEL1_DS:S1A_IW_GRDH_1SDV_20170907T103020_20170907...
    product:           GRDH
    safe:              S1A_IW_GRDH_1SDV_20170907T103020_20170907T103045_01826...
    swath:             IW
    multidataset:      False
    ...                ...
    start_date:        2017-09-07 10:30:20.936409
    stop_date:         2017-09-07 10:30:45.935264
    footprint:         POLYGON ((-67.84221143971432 20.72564283093837, -70.22...
    coverage:          170km * 251km (line * sample )
    orbit_pass:        Descending
    platform_heading:  -167.7668824808032

RadarSat2 example

[12]:

# get test file. You can replace with an path to other SAFE
filename = xsar.get_test_file('RS2_OK135107_PK1187782_DK1151894_SCWA_20220407_182127_VV_VH_SGF')

Open a dataset with a xsar.RadarSat2Meta object

A xsar.RadarSat2Meta object handles all attributes and methods that can’t be embdeded in a xarray.Dataset object. It can also replace a filename in xarray.open_dataset

[13]:

sar_meta = xsar.RadarSat2Meta(filename)
sar_meta

/home/datawork-cersat-public/cache/project/mpc-sentinel1/workspace/mamba/envs/testgeoviewsnumpy/lib/python3.11/site-packages/xradarsat2/radarSat2_xarray_reader.py:299: UserWarning: no explicit representation of timezones available for np.datetime64
  timestamp.append(np.datetime64(value["timeStamp"]).astype("datetime64[ns]"))
/home/datawork-cersat-public/cache/project/mpc-sentinel1/workspace/mamba/envs/testgeoviewsnumpy/lib/python3.11/site-packages/xradarsat2/radarSat2_xarray_reader.py:462: UserWarning: no explicit representation of timezones available for np.datetime64
  timestamp.append(np.datetime64(value["timeStamp"]).astype("datetime64[ns]"))
/home/datawork-cersat-public/cache/project/mpc-sentinel1/workspace/mamba/envs/testgeoviewsnumpy/lib/python3.11/site-packages/xradarsat2/radarSat2_xarray_reader.py:570: UserWarning: no explicit representation of timezones available for np.datetime64
  times.append(np.datetime64(value["timeOfDopplerCentroidEstimate"]).astype("datetime64[ns]"))
/home/datawork-cersat-public/cache/project/mpc-sentinel1/workspace/mamba/envs/testgeoviewsnumpy/lib/python3.11/site-packages/xradarsat2/radarSat2_xarray_reader.py:2056: UserWarning: no explicit representation of timezones available for np.datetime64
  final_dic[key] = np.datetime64(dic[key]).astype("datetime64[ns]")

[13]:

<RadarSat2Meta single object>

If holoviews extension is loaded, the <RadarSat2Meta objet> have a nice representation. (matplolib is also a valid extension)

[14]:

sar_meta

[14]:

<RadarSat2Meta single object>

sar_meta object is an xsar.RadarSat2Meta object that can be given to xarray.open_dataset or xsar.RadarSat2Dataset , as if it was a filename:

[15]:

sar_ds = xsar.RadarSat2Dataset(sar_meta)
sar_ds

/home/datawork-cersat-public/cache/project/mpc-sentinel1/workspace/mamba/envs/testgeoviewsnumpy/lib/python3.11/site-packages/xradarsat2/radarSat2_xarray_reader.py:299: UserWarning: no explicit representation of timezones available for np.datetime64
  timestamp.append(np.datetime64(value["timeStamp"]).astype("datetime64[ns]"))
/home/datawork-cersat-public/cache/project/mpc-sentinel1/workspace/mamba/envs/testgeoviewsnumpy/lib/python3.11/site-packages/xradarsat2/radarSat2_xarray_reader.py:462: UserWarning: no explicit representation of timezones available for np.datetime64
  timestamp.append(np.datetime64(value["timeStamp"]).astype("datetime64[ns]"))
/home/datawork-cersat-public/cache/project/mpc-sentinel1/workspace/mamba/envs/testgeoviewsnumpy/lib/python3.11/site-packages/xradarsat2/radarSat2_xarray_reader.py:570: UserWarning: no explicit representation of timezones available for np.datetime64
  times.append(np.datetime64(value["timeOfDopplerCentroidEstimate"]).astype("datetime64[ns]"))
/home/datawork-cersat-public/cache/project/mpc-sentinel1/workspace/mamba/envs/testgeoviewsnumpy/lib/python3.11/site-packages/xradarsat2/radarSat2_xarray_reader.py:2056: UserWarning: no explicit representation of timezones available for np.datetime64
  final_dic[key] = np.datetime64(dic[key]).astype("datetime64[ns]")

[15]:

<RadarSat2Dataset full coverage object>

Open a dataset at lower resolution

resolution keyword can be used to open a dataset at lower resolution.

It might be:

a dict {'line': 20, 'sample': 20}: 20*20 pixels. so if sensor resolution is 10m, the final resolution will be 100m
a string like '200m': Sensor resolution will be automatically used to compute the window size

Then we can instantiate a xsar.RadarSat2Dataset, with the given resolution. Note that the above pixel size has changed.

[16]:

sar_ds = xsar.RadarSat2Dataset(sar_meta, resolution='200m')
sar_ds

/home/datawork-cersat-public/cache/project/mpc-sentinel1/workspace/mamba/envs/testgeoviewsnumpy/lib/python3.11/site-packages/xradarsat2/radarSat2_xarray_reader.py:299: UserWarning: no explicit representation of timezones available for np.datetime64
  timestamp.append(np.datetime64(value["timeStamp"]).astype("datetime64[ns]"))
/home/datawork-cersat-public/cache/project/mpc-sentinel1/workspace/mamba/envs/testgeoviewsnumpy/lib/python3.11/site-packages/xradarsat2/radarSat2_xarray_reader.py:462: UserWarning: no explicit representation of timezones available for np.datetime64
  timestamp.append(np.datetime64(value["timeStamp"]).astype("datetime64[ns]"))
/home/datawork-cersat-public/cache/project/mpc-sentinel1/workspace/mamba/envs/testgeoviewsnumpy/lib/python3.11/site-packages/xradarsat2/radarSat2_xarray_reader.py:570: UserWarning: no explicit representation of timezones available for np.datetime64
  times.append(np.datetime64(value["timeOfDopplerCentroidEstimate"]).astype("datetime64[ns]"))
/home/datawork-cersat-public/cache/project/mpc-sentinel1/workspace/mamba/envs/testgeoviewsnumpy/lib/python3.11/site-packages/xradarsat2/radarSat2_xarray_reader.py:2056: UserWarning: no explicit representation of timezones available for np.datetime64
  final_dic[key] = np.datetime64(dic[key]).astype("datetime64[ns]")

[16]:

<RadarSat2Dataset full coverage object>

Extract a sub image of 10*10km around a lon/lat point

Convert (lon,lat) to (line, sample)

we can use sar_meta.ll2coords to convert (lon,lat) to (line, sample):

[17]:

# from a shapely object
point_lonlat =  sar_meta.footprint.centroid
point_coords = sar_meta.ll2coords(point_lonlat.x, point_lonlat.y)
point_coords

[17]:

(np.float64(5120.3748839621185), np.float64(5316.257595315285))

The result is floating, because it is the position inside the pixel. If real indexes from existing dataset is needed, you’ll have to use sar_ds.ll2coords Result will be the nearest (line, sample) in the dataset

[18]:

point_coords = sar_ds.ll2coords(point_lonlat.x, point_lonlat.y)
point_coords

[18]:

(5121.5, 5317.5)

Extract the sub-image

[19]:

box_size = 10000 # 10km
dist = {'line' : int(np.round(box_size / 2 / sar_meta.pixel_line_m)), 'sample': int(np.round(box_size / 2 / sar_meta.pixel_sample_m))}
dist

[19]:

{'line': 100, 'sample': 100}

The xarray/dask dataset is available as a property : sar_ds.dataset. This attribute can be set to a new values, so the attributes like pixel spacing and coverage are correctly recomputed:

[20]:

# select 10*10 km around point_coords
sar_ds.dataset = sar_ds.dataset.sel(line=slice(point_coords[0] - dist['line'], point_coords[0] + dist['line']), sample=slice(point_coords[1] - dist['sample'], point_coords[1] + dist['sample']))
sar_ds

[20]:

<RadarSat2Dataset sliced object>

[21]:

sar_ds.dataset

[21]:

<xarray.Dataset> Size: 1GB
Dimensions:          (line: 2569, pol: 2, sample: 2654)
Coordinates:
  * line             (line) float64 21kB 1.5 5.5 9.5 ... 1.027e+04 1.027e+04
  * pol              (pol) <U2 16B 'VV' 'VH'
  * sample           (sample) float64 21kB 1.5 5.5 9.5 ... 1.061e+04 1.061e+04
Data variables: (12/23)
    digital_number   (pol, line, sample) uint16 27MB dask.array<chunksize=(1, 2569, 2654), meta=np.ndarray>
    lines_flipped    bool 1B False
    samples_flipped  bool 1B True
    sampleSpacing    float64 8B 200.0
    lineSpacing      float64 8B 200.0
    sigma0_raw       (pol, line, sample) float64 109MB dask.array<chunksize=(1, 2569, 2654), meta=np.ndarray>
    ...               ...
    time             (line) datetime64[ns] 21kB dask.array<chunksize=(2569,), meta=np.ndarray>
    latitude         (line, sample) float64 55MB dask.array<chunksize=(2569, 2654), meta=np.ndarray>
    longitude        (line, sample) float64 55MB dask.array<chunksize=(2569, 2654), meta=np.ndarray>
    altitude         (line, sample) float64 55MB dask.array<chunksize=(2569, 2654), meta=np.ndarray>
    incidence        (line, sample) float64 55MB dask.array<chunksize=(2569, 2654), meta=np.ndarray>
    elevation        (line, sample) float64 55MB dask.array<chunksize=(2569, 2654), meta=np.ndarray>
Attributes: (12/18)
    product_path:           /home1/scratch/agrouaze/xsardatasync/xsardata/RS2...
    satellite:              RADARSAT-2
    inputDatasetId:         /Fred/RSAT-2/610044P
    rawDataStartTime:       2022-04-07T18:21:27.688416000
    satelliteHeight_units:  m
    satelliteHeight:        800612.0083192665
    ...                     ...
    stop_date:              2022-04-07 18:22:44.030964
    footprint:              POLYGON ((168.8147284622835 -19.82390747944243, 1...
    coverage:               515km * 530km (line * sample )
    pixel_line_m:           50.0
    pixel_sample_m:         50.0
    approx_transform:       |-0.00,-0.00, 168.85|\n|-0.00, 0.00,-19.85|\n| 0....

RCM example

[22]:

import warnings
warnings.filterwarnings('ignore')

[23]:

# get test file. You can replace with an path to other SAFE
filename = xsar.get_test_file('RCM1_OK1050603_PK1050605_1_SC50MB_20200214_115905_HH_HV_Z010')

Open a dataset with a xsar.RcmMeta object

A xsar.RcmMeta object handles all attributes and methods that can’t be embdeded in a xarray.Dataset object. It can also replace a filename in xarray.open_dataset

[24]:

sar_meta = xsar.RcmMeta(filename)
sar_meta

[24]:

<RcmMeta single object>

If holoviews extension is loaded, the <RcmMeta objet> have a nice representation. (matplolib is also a valid extension)

[25]:

sar_meta

[25]:

<RcmMeta single object>

sar_meta object is an xsar.RcmMeta object that can be given to xarray.open_dataset or xsar.RcmDataset , as if it was a filename:

[26]:

sar_ds = xsar.RcmDataset(sar_meta)
sar_ds

[26]:

<RcmDataset full coverage object>

Open a dataset at lower resolution

resolution keyword can be used to open a dataset at lower resolution.

It might be:

a dict {'line': 20, 'sample': 20}: 20*20 pixels. so if sensor resolution is 10m, the final resolution will be 100m
a string like '200m': Sensor resolution will be automatically used to compute the window size

Then we can instantiate a xsar.RcmDataset, with the given resolution. Note that the above pixel size has changed.

[27]:

sar_ds = xsar.RcmDataset(sar_meta, resolution='200m')
sar_ds

[27]:

<RcmDataset full coverage object>

Extract a sub image of 10*10km around a lon/lat point

Convert (lon,lat) to (line, sample)

we can use sar_meta.ll2coords to convert (lon,lat) to (line, sample):

[28]:

# from a shapely object
point_lonlat =  sar_meta.footprint.centroid
point_coords = sar_meta.ll2coords(point_lonlat.x, point_lonlat.y)
point_coords

[28]:

(np.float64(9601.981149464264), np.float64(9000.263581544685))

The result is floating, because it is the position inside the pixel. If real indexes from existing dataset is needed, you’ll have to use sar_ds.ll2coords Result will be the nearest (line, sample) in the dataset

[29]:

point_coords = sar_ds.ll2coords(point_lonlat.x, point_lonlat.y)
point_coords

[29]:

(9604.5, 9004.5)

Extract the sub-image

[30]:

box_size = 10000 # 10km
dist = {'line' : int(np.round(box_size / 2 / sar_meta.pixel_line_m)), 'sample': int(np.round(box_size / 2 / sar_meta.pixel_sample_m))}
dist

[30]:

{'line': 250, 'sample': 250}

The xarray/dask dataset is available as a property : sar_ds.dataset. This attribute can be set to a new values, so the attributes like pixel spacing and coverage are correctly recomputed:

[31]:

# select 10*10 km around point_coords
sar_ds.dataset = sar_ds.dataset.sel(line=slice(point_coords[0] - dist['line'], point_coords[0] + dist['line']), sample=slice(point_coords[1] - dist['sample'], point_coords[1] + dist['sample']))
sar_ds

[31]:

<RcmDataset sliced object>

[32]:

sar_ds.dataset

[32]:

<xarray.Dataset> Size: 569MB
Dimensions:          (line: 1912, pol: 2, sample: 1804)
Coordinates:
  * line             (line) float64 15kB 4.5 14.5 24.5 ... 1.91e+04 1.911e+04
  * pol              (pol) <U2 16B 'HH' 'HV'
  * sample           (sample) float64 14kB 4.5 14.5 24.5 ... 1.802e+04 1.803e+04
Data variables: (12/21)
    digital_number   (pol, line, sample) uint16 14MB dask.array<chunksize=(1, 1912, 1804), meta=np.ndarray>
    lines_flipped    bool 1B False
    samples_flipped  bool 1B True
    sampleSpacing    float64 8B 200.0
    lineSpacing      float64 8B 200.0
    sigma0_raw       (pol, line, sample) float64 55MB dask.array<chunksize=(1, 1912, 1804), meta=np.ndarray>
    ...               ...
    ground_heading   (line, sample) float32 14MB dask.array<chunksize=(1912, 1804), meta=np.ndarray>
    latitude         (line, sample) float64 28MB dask.array<chunksize=(1912, 1804), meta=np.ndarray>
    longitude        (line, sample) float64 28MB dask.array<chunksize=(1912, 1804), meta=np.ndarray>
    altitude         (line, sample) float64 28MB dask.array<chunksize=(1912, 1804), meta=np.ndarray>
    incidence        (line, sample) float32 14MB dask.array<chunksize=(1912, 1804), meta=np.ndarray>
    elevation        (line, sample) float64 28MB dask.array<chunksize=(1912, 1804), meta=np.ndarray>
Attributes: (12/18)
    name:                    RCM_DS:/home1/scratch/agrouaze/xsardatasync/xsar...
    short_name:              RCM_DS:RCM1_OK1050603_PK1050605_1_SC50MB_2020021...
    product:                 Z010
    safe:                    RCM1_OK1050603_PK1050605_1_SC50MB_20200214_11590...
    swath:                   SC50MB
    multidataset:            False
    ...                      ...
    stop_date:               2020-02-14 12:00:02.000000
    footprint:               POLYGON ((-84.63397773377159 48.71534850877886, ...
    coverage:                381km * 359km (line * sample )
    pixel_line_m:            20.0
    pixel_sample_m:          20.0
    approx_transform:        |-0.00,-0.00,-84.71|\n|-0.00, 0.00, 48.73|\n| 0....