Ctrl+K

RTC dataset comparison

RTC dataset comparison#

Now that we have constructed and organized data cubes for both the ASF and Planetary Computer datasets, it would be useful and informative to compare the two. This notebook will demonstrate steps to ensure a direct comparison between the two datasets as well as use of xarray plotting tools for a preliminary visual comparison.

Learning goals#

Xarray and python techniques:

  • Conditional selection based on non-dimensional coordinates using xr.Dataset.where()

  • Subsetting datasets based on dimensional coordinates using xr.DataArray.isin()

  • Adding dimensional and non-dimensional coordinates to xr.Dataset objects

  • Xarray plotting methods

  • Projecting xarray objects to different grids using xr.interp_like()

High-level science goals

  • Comparing and evaluating multiple datasets

  • Setting up multiple datasets for direct comparisons

  • Handling differences in spatial resolution

Other useful resources#

These are resources that contain additional examples and discussion of the content in this notebook and more.

%xmode minimal

import xarray as xr

import os
import xarray as xr
import rioxarray as rio
import geopandas as gpd
from datetime import datetime
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import dask

%matplotlib inline

Overview#

While the datasets we will be comparing in this notebook are quite similar, they contain important differences related to how they were processed. The product description pages for both datasets contain important information for understanding how they are generated.

Source data#

One important difference is the source data used to generate the RTC images. ASF uses the Single Look Complex (SLC) images that contains both amplitude and phase information for each pixel, is in radar coordinates and has not yet been multi-looked. The Microsoft Planetary Computer RTC imagery is generate from ground range detected images (GRD). GRD data has been detected, multi-looked and projected to ground range.

DEM#

A digital elevation model (DEM) is an important input parameter for RTC processing. ASF RTC processing uses the Copernicus DEM with 30 m resolution. Planetary Computer uses higher resolution DEMs. One effect of the use of different DEMs during RTC processing is the different spatial resolutions of the RTC products: Planetary Computer RTC imagery has a spatial resolution of 10 meters while ASF RTC imagery has a spatial resolution of 30 meters.

Software and setup#

import xarray as xr

import s1_tools

from s1_tools import points2coords

import os
import xarray as xr
import rioxarray as rio
import geopandas as gpd
from datetime import datetime
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import dask
import holoviews as hv
from holoviews import opts

Utility functions#

def power_to_db(input_arr):
    return (10*np.log10(np.abs(input_arr)))

def asf_pc_sidebyside(asf_input, pc_input, timestep):
    fig, axs = plt.subplots(ncols=2, figsize=(15,10))

    power_to_db(asf_input.vv.isel(acq_date=timestep)).plot(ax=axs[0], cmap=plt.cm.Greys_r, label = 'ASF');
    power_to_db(pc_input.sel(band='vv').isel(time=timestep)).plot(ax=axs[1], cmap=plt.cm.Greys_r, label = 'PC');

def single_time_mean_compare(asf_input, pc_input, time):
    fig, ax = plt.subplots(figsize=(8,8))
    power_to_db(asf_input['vv_'].isel(acq_date=time).mean(dim=['x','y'])).plot(ax=ax)
    power_to_db(pc_input.sel(band='vv').isel(time=time).mean(dim=['x','y'])).plot(ax=ax, color='red');

Read in prepared data#

We can use the storemagic command %store to retrieve the variable we constructed and saved in a previous notebook, rather than having to create it again. Read more about this here

This let’s us call the ASF dataset (vrt_new) and PC dataset (da_pc)

%store -r vrt_new da_pc

Read in the ASF dataset, which we named vrt_full in the earlier notebook:

vrt_new

<xarray.Dataset>
Dimensions:                           (acq_date: 95, y: 396, x: 290)
Coordinates: (12/20)
  * acq_date                          (acq_date) datetime64[ns] 2021-05-02 .....
    granule_id                        (acq_date) <U67 'S1A_IW_SLC__1SDV_20210...
  * x                                 (x) float64 6.194e+05 ... 6.281e+05
  * y                                 (y) float64 3.102e+06 ... 3.09e+06
    spatial_ref                       int64 0
    sensor                            (acq_date) <U3 'S1A' 'S1A' ... 'S1A' 'S1A'
    ...                                ...
    filtered                          (acq_date) <U1 'n' 'n' 'n' ... 'n' 'n' 'n'
    area                              (acq_date) <U1 'e' 'e' 'e' ... 'e' 'e' 'e'
    product_id                        (acq_date) <U4 '748F' '0D1E' ... 'BD36'
    acq_hour                          (acq_date) int64 0 12 12 12 12 ... 0 0 0 0
    orbital_dir                       (acq_date) <U4 'asc' 'desc' ... 'asc'
    data_take_id                      (acq_date) <U6 '047321' ... '052C00'
Data variables:
    vv                                (acq_date, y, x) float32 dask.array<chunksize=(1, 396, 290), meta=np.ndarray>
    vh                                (acq_date, y, x) float32 dask.array<chunksize=(1, 396, 290), meta=np.ndarray>
    ls                                (acq_date, y, x) float64 dask.array<chunksize=(1, 396, 290), meta=np.ndarray>

Before we go further, the nodata value for the ASF dataset is currently zero. Change all zero values to NaN:

vrt_new =  vrt_new.where(vrt_new.vv != 0., np.nan, drop=False)

da_pc

<xarray.DataArray 'stackstac-4d29467097b7d384d1232c9847300ccd' (time: 100,
                                                                band: 2,
                                                                y: 1188, x: 869)>
dask.array<fetch_raster_window, shape=(100, 2, 1188, 869), dtype=float64, chunksize=(1, 1, 1024, 869), chunktype=numpy.ndarray>
Coordinates: (12/40)
  * time                                   (time) datetime64[ns] 2021-06-02T1...
    id                                     (time) <U66 'S1A_IW_GRDH_1SDV_2021...
  * band                                   (band) <U2 'vh' 'vv'
  * x                                      (x) float64 6.194e+05 ... 6.281e+05
  * y                                      (y) float64 3.102e+06 ... 3.09e+06
    sar:looks_azimuth                      int64 1
    ...                                     ...
    s1:processing_level                    <U1 '1'
    title                                  (band) <U41 'VH: vertical transmit...
    description                            (band) <U173 'Terrain-corrected ga...
    raster:bands                           object {'nodata': -32768, 'data_ty...
    epsg                                   int64 32645
    granule_id                             (time) <U62 'S1A_IW_GRDH_1SDV_2021...
Attributes:
    spec:        RasterSpec(epsg=32645, bounds=(619420.0, 3089780.0, 628110.0...
    crs:         epsg:32645
    transform:   | 10.00, 0.00, 619420.00|\n| 0.00,-10.00, 3101660.00|\n| 0.0...
    resolution:  10.0

Extract common data take ID from granule IDs#

We want to ensure that we are performing a direct comparison of the ASF and PC datasets. To do this, we would like to use the acquisition ID that is stored in the source granule name (published by ESA). In the setup notebooks we attached the entire granule IDs of the SLC images to the ASF dataset and the GRD images to the PC dataset. In the ASF data inspection notebook, we attached data take id as a non-dimensional coordinate. Now we will do the same for the Planetary Computer dataset, extracting just the 6-digit acquisition ID from the granule ID and using this for a scene-by-scene comparison.

data_take_pc = [str(da_pc.isel(time=t).granule_id.values)[56:] for t in range(len(da_pc.time))]

Assign data_take_id as a non-dimensional coordinate. Rather than use the xr.assign_coords() function, we are adding the object by specifying a tuple with the form (‘coord_name’, coord_data, attrs):

da_pc.coords['data_take_id'] = ('time', data_take_pc)

vrt_new

<xarray.Dataset>
Dimensions:                           (acq_date: 95, y: 396, x: 290)
Coordinates: (12/20)
  * acq_date                          (acq_date) datetime64[ns] 2021-05-02 .....
    granule_id                        (acq_date) <U67 'S1A_IW_SLC__1SDV_20210...
  * x                                 (x) float64 6.194e+05 ... 6.281e+05
  * y                                 (y) float64 3.102e+06 ... 3.09e+06
    spatial_ref                       int64 0
    sensor                            (acq_date) <U3 'S1A' 'S1A' ... 'S1A' 'S1A'
    ...                                ...
    filtered                          (acq_date) <U1 'n' 'n' 'n' ... 'n' 'n' 'n'
    area                              (acq_date) <U1 'e' 'e' 'e' ... 'e' 'e' 'e'
    product_id                        (acq_date) <U4 '748F' '0D1E' ... 'BD36'
    acq_hour                          (acq_date) int64 0 12 12 12 12 ... 0 0 0 0
    orbital_dir                       (acq_date) <U4 'asc' 'desc' ... 'asc'
    data_take_id                      (acq_date) <U6 '047321' ... '052C00'
Data variables:
    vv                                (acq_date, y, x) float32 dask.array<chunksize=(1, 396, 290), meta=np.ndarray>
    vh                                (acq_date, y, x) float32 dask.array<chunksize=(1, 396, 290), meta=np.ndarray>
    ls                                (acq_date, y, x) float64 dask.array<chunksize=(1, 396, 290), meta=np.ndarray>
da_pc

<xarray.DataArray 'stackstac-4d29467097b7d384d1232c9847300ccd' (time: 100,
                                                                band: 2,
                                                                y: 1188, x: 869)>
dask.array<fetch_raster_window, shape=(100, 2, 1188, 869), dtype=float64, chunksize=(1, 1, 1024, 869), chunktype=numpy.ndarray>
Coordinates: (12/41)
  * time                                   (time) datetime64[ns] 2021-06-02T1...
    id                                     (time) <U66 'S1A_IW_GRDH_1SDV_2021...
  * band                                   (band) <U2 'vh' 'vv'
  * x                                      (x) float64 6.194e+05 ... 6.281e+05
  * y                                      (y) float64 3.102e+06 ... 3.09e+06
    sar:looks_azimuth                      int64 1
    ...                                     ...
    title                                  (band) <U41 'VH: vertical transmit...
    description                            (band) <U173 'Terrain-corrected ga...
    raster:bands                           object {'nodata': -32768, 'data_ty...
    epsg                                   int64 32645
    granule_id                             (time) <U62 'S1A_IW_GRDH_1SDV_2021...
    data_take_id                           (time) <U6 '0480FD' ... '054B39'
Attributes:
    spec:        RasterSpec(epsg=32645, bounds=(619420.0, 3089780.0, 628110.0...
    crs:         epsg:32645
    transform:   | 10.00, 0.00, 619420.00|\n| 0.00,-10.00, 3101660.00|\n| 0.0...
    resolution:  10.0

Subset to only common data takes#

Now we have data take ID coordinates for both datasets. We want to find the common data take IDs between the two datasets. To do this, I extract a list of the acquisition IDs (data_take_id) for both datasets and then find the intersection of the two lists (the list object common_data_takes)

pc_data_take_ls = list(da_pc.data_take_id.values)
asf_data_take_ls = list(vrt_new.data_take_id.values)

common_data_takes  = list(set(pc_data_take_ls) & set(asf_data_take_ls))
len(common_data_takes)

83

It looks like there are 83 RTC images that are generated from common acquisitions between the two datasets

We’d like to subset the ASF and the PC datasets to only the common acquisitions. The xarray .isin() method is very useful for this type of selection.

subset_condition_asf = vrt_new.data_take_id.isin(common_data_takes)

subset_condition_pc = da_pc.data_take_id.isin(common_data_takes)

Next, we want to select the elements along the acq_ate (and time) dimensions where data_take_id satisfies the requirement of being in the common_data_takes list. I found this stack overflow answer useful for this step: https://stackoverflow.com/questions/70777676/xarray-select-dataarray-according-to-an-non-dimension-coordinate.

asf_subset = vrt_new.sel(acq_date = vrt_new.data_take_id.isin(common_data_takes))

pc_subset = da_pc.sel(time = da_pc.data_take_id.isin(common_data_takes))

pc_subset

<xarray.DataArray 'stackstac-4d29467097b7d384d1232c9847300ccd' (time: 83,
                                                                band: 2,
                                                                y: 1188, x: 869)>
dask.array<getitem, shape=(83, 2, 1188, 869), dtype=float64, chunksize=(1, 1, 1024, 869), chunktype=numpy.ndarray>
Coordinates: (12/41)
  * time                                   (time) datetime64[ns] 2021-06-02T1...
    id                                     (time) <U66 'S1A_IW_GRDH_1SDV_2021...
  * band                                   (band) <U2 'vh' 'vv'
  * x                                      (x) float64 6.194e+05 ... 6.281e+05
  * y                                      (y) float64 3.102e+06 ... 3.09e+06
    sar:looks_azimuth                      int64 1
    ...                                     ...
    title                                  (band) <U41 'VH: vertical transmit...
    description                            (band) <U173 'Terrain-corrected ga...
    raster:bands                           object {'nodata': -32768, 'data_ty...
    epsg                                   int64 32645
    granule_id                             (time) <U62 'S1A_IW_GRDH_1SDV_2021...
    data_take_id                           (time) <U6 '0480FD' ... '052C00'
Attributes:
    spec:        RasterSpec(epsg=32645, bounds=(619420.0, 3089780.0, 628110.0...
    crs:         epsg:32645
    transform:   | 10.00, 0.00, 619420.00|\n| 0.00,-10.00, 3101660.00|\n| 0.0...
    resolution:  10.0
asf_subset

<xarray.Dataset>
Dimensions:                           (acq_date: 83, y: 396, x: 290)
Coordinates: (12/20)
  * acq_date                          (acq_date) datetime64[ns] 2021-06-02 .....
    granule_id                        (acq_date) <U67 'S1A_IW_SLC__1SDV_20210...
  * x                                 (x) float64 6.194e+05 ... 6.281e+05
  * y                                 (y) float64 3.102e+06 ... 3.09e+06
    spatial_ref                       int64 0
    sensor                            (acq_date) <U3 'S1A' 'S1A' ... 'S1A' 'S1A'
    ...                                ...
    filtered                          (acq_date) <U1 'n' 'n' 'n' ... 'n' 'n' 'n'
    area                              (acq_date) <U1 'e' 'e' 'e' ... 'e' 'e' 'e'
    product_id                        (acq_date) <U4 '1424' 'ABA0' ... 'BD36'
    acq_hour                          (acq_date) int64 12 12 0 12 0 ... 0 0 0 0
    orbital_dir                       (acq_date) <U4 'desc' 'desc' ... 'asc'
    data_take_id                      (acq_date) <U6 '0480FD' ... '052C00'
Data variables:
    vv                                (acq_date, y, x) float32 dask.array<chunksize=(1, 396, 290), meta=np.ndarray>
    vh                                (acq_date, y, x) float32 dask.array<chunksize=(1, 396, 290), meta=np.ndarray>
    ls                                (acq_date, y, x) float64 dask.array<chunksize=(1, 396, 290), meta=np.ndarray>

Visual comparison#

Let’s perform some preliminary visual comparisons of the datasets:

Individual time steps#

Use asf_bc_sidebyside() to compare snapshot views of backscatter at individual time steps:

asf_pc_sidebyside(asf_subset, pc_subset, 2)
../_images/2592b8ceaf6c1a594034036e1817440d97ede543c444a1098fb82b779acb3dfa.png 
asf_pc_sidebyside(asf_subset, pc_subset, 4)
../_images/7414e079129739c5356e6a19f5e4f65f1f85b0f2556bcb56fd55c9257f69b931.png 
asf_pc_sidebyside(asf_subset, pc_subset, 7)
../_images/e9de4af3765472119f72bb8e243f6a366f12b818dabd29818152bf87787ee0da.png 
asf_pc_sidebyside(asf_subset, pc_subset, 12)
../_images/17b1813e78366a6ac0c437ef503894f9b4978c5c02348bd4f9a8d4b2d008e2f1.png

Average backscatter over time#

Now we will look at the backscatter mean along the time dimension

fig, ax = plt.subplots(figsize=(12,8))

power_to_db(asf_subset.vv.mean(dim=['x','y'])).plot(ax=ax, color='blue', linestyle='None', marker='o', markersize=5, label = 'ASF')
power_to_db(pc_subset.sel(band='vv').mean(dim=['x','y'])).plot(ax=ax, color='red', linestyle='None', marker='o', markersize=5, label='PC');

[<matplotlib.lines.Line2D>]
../_images/e85328617ec09a4085cb3c8d4f42618c1202d4dac89a0d10d552379d60683bf3.png

Handling spatial resolution differences#

The PC dataset has a higher spatial resolution than the ASF dataset, which you can see in the sizes of the x and y dimensions of the dataset. The ASF dataset masks out more pixels due to shadow than the PC dataset. To explore whether or not this is the cause of the observed backscatter offset between the two, we need to mask out the pixels in the PC dataset that are masked in the ASF dataset.

In order to do this, first we must downsample the PC dataset to match the spatial resolution of the ASF dataset. There are multiple xarray methods that could be used but here we will use the xr.interp_like() method because we would like the PC x and y dimensions to match the ASF dimensiosn exactly.

pc_downsample = pc_subset.interp_like(asf_subset)

pc_downsample

<xarray.DataArray 'stackstac-4d29467097b7d384d1232c9847300ccd' (time: 83,
                                                                band: 2,
                                                                y: 396, x: 290)>
dask.array<transpose, shape=(83, 2, 396, 290), dtype=float64, chunksize=(1, 1, 396, 290), chunktype=numpy.ndarray>
Coordinates: (12/41)
  * time                                   (time) datetime64[ns] 2021-06-02T1...
    id                                     (time) <U66 'S1A_IW_GRDH_1SDV_2021...
  * band                                   (band) <U2 'vh' 'vv'
    sar:looks_azimuth                      int64 1
    sat:relative_orbit                     (time) int64 114 12 48 ... 12 114 12
    s1:instrument_configuration_ID         (time) <U1 '6' '6' '6' ... '7' '7'
    ...                                     ...
    raster:bands                           object {'nodata': -32768, 'data_ty...
    epsg                                   int64 32645
    granule_id                             (time) <U62 'S1A_IW_GRDH_1SDV_2021...
    data_take_id                           (time) <U6 '0480FD' ... '052C00'
  * x                                      (x) float64 6.194e+05 ... 6.281e+05
  * y                                      (y) float64 3.102e+06 ... 3.09e+06
Attributes:
    spec:        RasterSpec(epsg=32645, bounds=(619420.0, 3089780.0, 628110.0...
    crs:         epsg:32645
    transform:   | 10.00, 0.00, 619420.00|\n| 0.00,-10.00, 3101660.00|\n| 0.0...
    resolution:  10.0
pc_mask = xr.where(asf_subset.notnull(), pc_downsample,np.nan)

/home/emmamarshall/miniconda3/envs/sentinel/lib/python3.10/site-packages/dask/array/core.py:4796: PerformanceWarning: Increasing number of chunks by factor of 83
  result = blockwise(
/home/emmamarshall/miniconda3/envs/sentinel/lib/python3.10/site-packages/dask/array/core.py:4796: PerformanceWarning: Increasing number of chunks by factor of 83
  result = blockwise(
/home/emmamarshall/miniconda3/envs/sentinel/lib/python3.10/site-packages/dask/array/core.py:4796: PerformanceWarning: Increasing number of chunks by factor of 83
  result = blockwise(

pc_mask

<xarray.Dataset>
Dimensions:                                (acq_date: 83, y: 396, x: 290,
                                            time: 83, band: 2)
Coordinates: (12/55)
  * acq_date                               (acq_date) datetime64[ns] 2021-06-...
  * x                                      (x) float64 6.194e+05 ... 6.281e+05
  * y                                      (y) float64 3.102e+06 ... 3.09e+06
    spatial_ref                            int64 0
    sensor                                 (acq_date) <U3 'S1A' 'S1A' ... 'S1A'
    beam_mode                              (acq_date) <U2 'IW' 'IW' ... 'IW'
    ...                                     ...
    s1:total_slices                        (time) <U2 '21' '20' ... '21' '20'
    s1:processing_level                    <U1 '1'
    title                                  (band) <U41 'VH: vertical transmit...
    description                            (band) <U173 'Terrain-corrected ga...
    raster:bands                           object {'nodata': -32768, 'data_ty...
    epsg                                   int64 32645
Data variables:
    vv                                     (acq_date, y, x, time, band) float64 dask.array<chunksize=(1, 396, 290, 1, 1), meta=np.ndarray>
    vh                                     (acq_date, y, x, time, band) float64 dask.array<chunksize=(1, 396, 290, 1, 1), meta=np.ndarray>
    ls                                     (acq_date, y, x, time, band) float64 dask.array<chunksize=(1, 396, 290, 1, 1), meta=np.ndarray>

Oops, this isn’t what we want. Now we have a 5-dimensional object instead of a 4-dimensional object, because the dimensions of the asf_subset object were broadcast onto the pc_downsample object, meaning that we have both a time and a acq_date object.

Let’s go back and try to fix this…

pc_downsample.isel(time=0).time.values

numpy.datetime64('2021-06-02T12:05:57.441074000')

asf_subset = asf_subset.rename({'acq_date':'time'})

asf_subset.isel(time=0).time.values

numpy.datetime64('2021-06-02T00:00:00.000000000')

Okay, so we need to take the time off of the pc_downsample time coordinates:

pc_downsample.coords['time'] = pc_downsample.time.dt.date

pc_downsample.coords.dtypes

Frozen({'time': dtype('O'), 'id': dtype('<U66'), 'band': dtype('<U2'), 'sar:looks_azimuth': dtype('int64'), 'sat:relative_orbit': dtype('int64'), 's1:instrument_configuration_ID': dtype('<U1'), 's1:orbit_source': dtype('<U8'), 'end_datetime': dtype('<U32'), 'sat:orbit_state': dtype('<U10'), 'platform': dtype('<U11'), 's1:slice_number': dtype('<U2'), 'start_datetime': dtype('<U32'), 'sar:resolution_azimuth': dtype('int64'), 's1:resolution': dtype('<U4'), 'sar:frequency_band': dtype('<U1'), 'sar:center_frequency': dtype('float64'), 'sar:pixel_spacing_azimuth': dtype('int64'), 'sat:platform_international_designator': dtype('<U9'), 'constellation': dtype('<U10'), 'sar:polarizations': dtype('O'), 'proj:epsg': dtype('int64'), 'sar:pixel_spacing_range': dtype('int64'), 'sar:looks_range': dtype('int64'), 'sar:resolution_range': dtype('int64'), 'sat:absolute_orbit': dtype('int64'), 'sar:looks_equivalent_number': dtype('float64'), 'sar:observation_direction': dtype('<U5'), 'sar:product_type': dtype('<U3'), 's1:product_timeliness': dtype('<U8'), 'sar:instrument_mode': dtype('<U2'), 's1:datatake_id': dtype('<U6'), 's1:total_slices': dtype('<U2'), 's1:processing_level': dtype('<U1'), 'title': dtype('<U41'), 'description': dtype('<U173'), 'raster:bands': dtype('O'), 'epsg': dtype('int64'), 'granule_id': dtype('<U62'), 'data_take_id': dtype('<U6'), 'x': dtype('float64'), 'y': dtype('float64')})

asf_subset.coords.dtypes

Frozen({'time': dtype('<M8[ns]'), 'granule_id': dtype('<U67'), 'x': dtype('float64'), 'y': dtype('float64'), 'spatial_ref': dtype('int64'), 'sensor': dtype('<U3'), 'beam_mode': dtype('<U2'), 'acquisition_time': dtype('<M8[ns]'), 'polarisation_type': dtype('<U2'), 'orbit_type': dtype('<U1'), 'terrain_correction_pixel_spacing': dtype('<U5'), 'output_format': dtype('<U1'), 'output_type': dtype('<U1'), 'masked': dtype('<U1'), 'filtered': dtype('<U1'), 'area': dtype('<U1'), 'product_id': dtype('<U4'), 'acq_hour': dtype('int64'), 'orbital_dir': dtype('<U4'), 'data_take_id': dtype('<U6')})

Now let’s try again:

pc_mask = xr.where(asf_subset.notnull(), pc_downsample,np.nan)

pc_mask

<xarray.Dataset>
Dimensions:                                (time: 83, y: 396, x: 290, band: 2)
Coordinates: (12/55)
  * time                                   (time) datetime64[ns] 2021-06-02 ....
  * x                                      (x) float64 6.194e+05 ... 6.281e+05
  * y                                      (y) float64 3.102e+06 ... 3.09e+06
    spatial_ref                            int64 0
    sensor                                 (time) <U3 'S1A' 'S1A' ... 'S1A'
    beam_mode                              (time) <U2 'IW' 'IW' ... 'IW' 'IW'
    ...                                     ...
    s1:total_slices                        (time) <U2 '21' '20' ... '21' '20'
    s1:processing_level                    <U1 '1'
    title                                  (band) <U41 'VH: vertical transmit...
    description                            (band) <U173 'Terrain-corrected ga...
    raster:bands                           object {'nodata': -32768, 'data_ty...
    epsg                                   int64 32645
Data variables:
    vv                                     (time, y, x, band) float64 dask.array<chunksize=(1, 396, 290, 1), meta=np.ndarray>
    vh                                     (time, y, x, band) float64 dask.array<chunksize=(1, 396, 290, 1), meta=np.ndarray>
    ls                                     (time, y, x, band) float64 dask.array<chunksize=(1, 396, 290, 1), meta=np.ndarray>

Great! Looks like its working

asf_subset

<xarray.Dataset>
Dimensions:                           (time: 83, y: 396, x: 290)
Coordinates: (12/20)
  * time                              (time) datetime64[ns] 2021-06-02 ... 20...
    granule_id                        (time) <U67 'S1A_IW_SLC__1SDV_20210602T...
  * x                                 (x) float64 6.194e+05 ... 6.281e+05
  * y                                 (y) float64 3.102e+06 ... 3.09e+06
    spatial_ref                       int64 0
    sensor                            (time) <U3 'S1A' 'S1A' ... 'S1A' 'S1A'
    ...                                ...
    filtered                          (time) <U1 'n' 'n' 'n' 'n' ... 'n' 'n' 'n'
    area                              (time) <U1 'e' 'e' 'e' 'e' ... 'e' 'e' 'e'
    product_id                        (time) <U4 '1424' 'ABA0' ... '5BAD' 'BD36'
    acq_hour                          (time) int64 12 12 0 12 0 ... 12 0 0 0 0
    orbital_dir                       (time) <U4 'desc' 'desc' ... 'asc' 'asc'
    data_take_id                      (time) <U6 '0480FD' '048318' ... '052C00'
Data variables:
    vv                                (time, y, x) float32 dask.array<chunksize=(1, 396, 290), meta=np.ndarray>
    vh                                (time, y, x) float32 dask.array<chunksize=(1, 396, 290), meta=np.ndarray>
    ls                                (time, y, x) float64 dask.array<chunksize=(1, 396, 290), meta=np.ndarray>
pc_mask.sel(band='vv').mean(dim=['x','y'])

<xarray.Dataset>
Dimensions:                                (time: 83)
Coordinates: (12/53)
  * time                                   (time) datetime64[ns] 2021-06-02 ....
    spatial_ref                            int64 0
    sensor                                 (time) <U3 'S1A' 'S1A' ... 'S1A'
    beam_mode                              (time) <U2 'IW' 'IW' ... 'IW' 'IW'
    acquisition_time                       (time) datetime64[ns] 2022-10-21T1...
    polarisation_type                      (time) <U2 'DV' 'DV' ... 'DV' 'DV'
    ...                                     ...
    s1:total_slices                        (time) <U2 '21' '20' ... '21' '20'
    s1:processing_level                    <U1 '1'
    title                                  <U41 'VV: vertical transmit, verti...
    description                            <U173 'Terrain-corrected gamma nau...
    raster:bands                           object {'nodata': -32768, 'data_ty...
    epsg                                   int64 32645
Data variables:
    vv                                     (time) float64 dask.array<chunksize=(1,), meta=np.ndarray>
    vh                                     (time) float64 dask.array<chunksize=(1,), meta=np.ndarray>
    ls                                     (time) float64 dask.array<chunksize=(1,), meta=np.ndarray>
fig, ax = plt.subplots(figsize=(12,8))

power_to_db(asf_subset.vv.mean(dim=['x','y'])).plot(ax=ax, color='blue', linestyle='None', marker='o', markersize=5, label = 'ASF');
power_to_db(pc_mask.sel(band='vv').mean(dim=['x','y'], skipna=True).vv).plot(x = 'time',ax=ax, color='red', linestyle='None', marker='o', markersize=5, label='PC');
../_images/dae808fa50e639f4b18b67220545ec5245316b3c979daef83f4600059038dbee.png

Break up by asc, desc pass#

Use xr.where() to subset datasets by orbital pass:

asf_subset

<xarray.Dataset>
Dimensions:                           (acq_date: 90, y: 396, x: 290)
Coordinates: (12/20)
  * acq_date                          (acq_date) datetime64[ns] 2021-06-02 .....
    granule_id                        (acq_date) <U67 'S1A_IW_SLC__1SDV_20210...
  * x                                 (x) float64 6.194e+05 ... 6.281e+05
  * y                                 (y) float64 3.102e+06 ... 3.09e+06
    spatial_ref                       int64 0
    sensor                            (acq_date) <U3 'S1A' 'S1A' ... 'S1A' 'S1A'
    ...                                ...
    filtered                          (acq_date) <U1 'n' 'n' 'n' ... 'n' 'n' 'n'
    area                              (acq_date) <U1 'e' 'e' 'e' ... 'e' 'e' 'e'
    product_id                        (acq_date) <U4 '1424' 'ABA0' ... 'BD36'
    acq_hour                          (acq_date) int64 12 12 0 12 ... 0 12 12 0
    orbital_dir                       (acq_date) <U4 'desc' 'desc' ... 'asc'
    data_take_id                      (acq_date) <U6 '0480FD' ... '052C00'
Data variables:
    vv                                (acq_date, y, x) float32 dask.array<chunksize=(1, 396, 290), meta=np.ndarray>
    vh                                (acq_date, y, x) float32 dask.array<chunksize=(1, 396, 290), meta=np.ndarray>
    ls                                (acq_date, y, x) float64 dask.array<chunksize=(1, 396, 290), meta=np.ndarray>
asf_desc = asf_subset.where(asf_subset['orbital_dir'] == 'desc', drop=True)
asf_asc = asf_subset.where(asf_subset['orbital_dir'] == 'asc', drop=True)

#asf_asc.vv.mean(dim=['x','y']).plot(marker='o', linestyle='None')

pc_desc = pc_subset.where(pc_subset['sat:orbit_state'] == 'descending', drop=True)
pc_asc = pc_subset.where(pc_subset['sat:orbit_state'] == 'ascending', drop=True)

fig, axs = plt.subplots(ncols =2 ,figsize=(16,8))

power_to_db(asf_desc.vv.mean(dim=['x','y'])).plot(ax=axs[0], color = 'blue', linestyle='None', marker='+', markersize=5, label = 'ASF VV desc');
power_to_db(pc_desc.sel(band='vv').mean(dim=['x','y'])).plot(ax=axs[0], color = 'red', linestyle='None', marker = '+', markersize=5, label = 'PC VV desc');
            
power_to_db(asf_asc.vv.mean(dim=['x','y'])).plot(ax=axs[0], color = 'blue', linestyle='None', marker='o', markersize=5, label = 'ASF VV asc');
power_to_db(pc_asc.sel(band='vv').mean(dim=['x','y'])).plot(ax=axs[0], color = 'red', linestyle='None', marker = 'o', markersize=5, label = 'PC VV asc');
axs[0].legend(loc = 'lower right')
axs[0].set_xlabel('Acquisition date')
axs[0].set_ylabel('Backscatter (sigma-nought) (dB)') 
axs[0].set_title('VV')

power_to_db(asf_desc.vh.mean(dim=['x','y'])).plot(ax=axs[1], color = 'blue', linestyle='None', marker='+', markersize=5, label = 'ASF VH desc');
power_to_db(pc_desc.sel(band='vh').mean(dim=['x','y'])).plot(ax=axs[1], color = 'red', linestyle='None', marker = '+', markersize=5, label = 'PC VH desc');
            
power_to_db(asf_asc.vh.mean(dim=['x','y'])).plot(ax=axs[1], color = 'blue', linestyle='None', marker='o', markersize=5, label = 'ASF VH asc');
power_to_db(pc_asc.sel(band='vh').mean(dim=['x','y'])).plot(ax=axs[1], color = 'red', linestyle='None', marker = 'o', markersize=5, label = 'PC VH asc');
axs[1].legend(loc = 'lower right')
axs[1].set_xlabel('Acquisition date')
axs[1].set_ylabel('Backscatter (sigma-nought) (dB)') 
axs[1].set_title('VH')

fig.suptitle('Comparison of backscatter values from ASF-processed and PC-processed Sentinel-1 RTC imagery');

Text(0.5, 0.98, 'Comparison of backscatter values from ASF-processed and PC-processed Sentinel-1 RTC imagery')
../_images/5eefdd52d326481b5f849595d58270bce9d079c33fc9cd353f7734122f5e577a.png
Contents