RTC time series analysis

application example: temporal variability in backscatter over Himalayan proglacial lakes

import geopandas as gpd
import xarray as xr
from s1_tools import power_to_db, get_bbox_single

Read in RTC data

• this example will use ASF dataset

%store -r

asf_clip

<xarray.Dataset>
Dimensions:      (x: 290, y: 396, acq_date: 103)
Coordinates:
  * x            (x) float64 6.194e+05 6.195e+05 ... 6.281e+05 6.281e+05
  * y            (y) float64 3.09e+06 3.09e+06 3.09e+06 ... 3.102e+06 3.102e+06
  * acq_date     (acq_date) datetime64[ns] 2021-05-02 2021-05-05 ... 2022-05-21
    band         int64 1
    spatial_ref  int64 0
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) float32 dask.array<chunksize=(1, 396, 290), meta=np.ndarray>

xarray.Dataset

• Dimensions:
  • x: 290
  • y: 396
  • acq_date: 103
• Coordinates: (5)
  • x
    (x)
    float64
    6.194e+05 6.195e+05 ... 6.281e+05

    axis :

    X

    long_name :

    x coordinate of projection

    standard_name :

    projection_x_coordinate

    units :

    metre

    array([619425., 619455., 619485., ..., 628035., 628065., 628095.])

  • y
    (y)
    float64
    3.09e+06 3.09e+06 ... 3.102e+06

    axis :

    Y

    long_name :

    y coordinate of projection

    standard_name :

    projection_y_coordinate

    units :

    metre

    array([3089805., 3089835., 3089865., ..., 3101595., 3101625., 3101655.])

  • acq_date
    (acq_date)
    datetime64[ns]
    2021-05-02 ... 2022-05-21

    array(['2021-05-02T00:00:00.000000000', '2021-05-05T00:00:00.000000000',
           '2021-05-09T00:00:00.000000000', '2021-05-14T00:00:00.000000000',
           '2021-05-14T00:00:00.000000000', '2021-05-17T00:00:00.000000000',
           '2021-05-21T00:00:00.000000000', '2021-05-26T00:00:00.000000000',
           '2021-05-29T00:00:00.000000000', '2021-06-02T00:00:00.000000000',
           '2021-06-07T00:00:00.000000000', '2021-06-10T00:00:00.000000000',
           '2021-06-14T00:00:00.000000000', '2021-06-19T00:00:00.000000000',
           '2021-06-22T00:00:00.000000000', '2021-06-26T00:00:00.000000000',
           '2021-07-01T00:00:00.000000000', '2021-07-01T00:00:00.000000000',
           '2021-07-04T00:00:00.000000000', '2021-07-08T00:00:00.000000000',
           '2021-07-13T00:00:00.000000000', '2021-07-13T00:00:00.000000000',
           '2021-07-13T00:00:00.000000000', '2021-07-16T00:00:00.000000000',
           '2021-07-20T00:00:00.000000000', '2021-07-25T00:00:00.000000000',
           '2021-07-28T00:00:00.000000000', '2021-08-01T00:00:00.000000000',
           '2021-08-06T00:00:00.000000000', '2021-08-09T00:00:00.000000000',
           '2021-08-13T00:00:00.000000000', '2021-08-18T00:00:00.000000000',
           '2021-08-21T00:00:00.000000000', '2021-08-25T00:00:00.000000000',
           '2021-08-30T00:00:00.000000000', '2021-09-02T00:00:00.000000000',
           '2021-09-06T00:00:00.000000000', '2021-09-11T00:00:00.000000000',
           '2021-09-14T00:00:00.000000000', '2021-09-18T00:00:00.000000000',
           '2021-09-23T00:00:00.000000000', '2021-09-26T00:00:00.000000000',
           '2021-09-30T00:00:00.000000000', '2021-10-05T00:00:00.000000000',
           '2021-10-08T00:00:00.000000000', '2021-10-12T00:00:00.000000000',
           '2021-10-17T00:00:00.000000000', '2021-10-20T00:00:00.000000000',
           '2021-10-24T00:00:00.000000000', '2021-10-29T00:00:00.000000000',
           '2021-11-01T00:00:00.000000000', '2021-11-05T00:00:00.000000000',
           '2021-11-10T00:00:00.000000000', '2021-11-13T00:00:00.000000000',
           '2021-11-17T00:00:00.000000000', '2021-11-22T00:00:00.000000000',
           '2021-11-25T00:00:00.000000000', '2021-11-29T00:00:00.000000000',
           '2021-12-04T00:00:00.000000000', '2021-12-07T00:00:00.000000000',
           '2021-12-11T00:00:00.000000000', '2021-12-16T00:00:00.000000000',
           '2021-12-19T00:00:00.000000000', '2021-12-23T00:00:00.000000000',
           '2021-12-28T00:00:00.000000000', '2021-12-31T00:00:00.000000000',
           '2022-01-04T00:00:00.000000000', '2022-01-09T00:00:00.000000000',
           '2022-01-12T00:00:00.000000000', '2022-01-16T00:00:00.000000000',
           '2022-01-21T00:00:00.000000000', '2022-01-24T00:00:00.000000000',
           '2022-01-28T00:00:00.000000000', '2022-02-02T00:00:00.000000000',
           '2022-02-05T00:00:00.000000000', '2022-02-09T00:00:00.000000000',
           '2022-02-14T00:00:00.000000000', '2022-02-17T00:00:00.000000000',
           '2022-02-21T00:00:00.000000000', '2022-02-26T00:00:00.000000000',
           '2022-03-01T00:00:00.000000000', '2022-03-05T00:00:00.000000000',
           '2022-03-10T00:00:00.000000000', '2022-03-13T00:00:00.000000000',
           '2022-03-17T00:00:00.000000000', '2022-03-22T00:00:00.000000000',
           '2022-03-25T00:00:00.000000000', '2022-03-29T00:00:00.000000000',
           '2022-04-03T00:00:00.000000000', '2022-04-03T00:00:00.000000000',
           '2022-04-06T00:00:00.000000000', '2022-04-10T00:00:00.000000000',
           '2022-04-15T00:00:00.000000000', '2022-04-18T00:00:00.000000000',
           '2022-04-22T00:00:00.000000000', '2022-04-27T00:00:00.000000000',
           '2022-04-30T00:00:00.000000000', '2022-05-04T00:00:00.000000000',
           '2022-05-09T00:00:00.000000000', '2022-05-16T00:00:00.000000000',
           '2022-05-21T00:00:00.000000000', '2022-05-21T00:00:00.000000000',
           '2022-05-21T00:00:00.000000000'], dtype='datetime64[ns]')

  • band
    ()
    int64
    1

    array(1)

  • spatial_ref
    ()
    int64
    0

    crs_wkt :

    PROJCS["WGS 84 / UTM zone 45N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",87],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32645"]]

    semi_major_axis :

    6378137.0

    semi_minor_axis :

    6356752.314245179

    inverse_flattening :

    298.257223563

    reference_ellipsoid_name :

    WGS 84

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    geographic_crs_name :

    WGS 84

    horizontal_datum_name :

    World Geodetic System 1984

    projected_crs_name :

    WGS 84 / UTM zone 45N

    grid_mapping_name :

    transverse_mercator

    latitude_of_projection_origin :

    0.0

    longitude_of_central_meridian :

    87.0

    false_easting :

    500000.0

    false_northing :

    0.0

    scale_factor_at_central_meridian :

    0.9996

    spatial_ref :

    PROJCS["WGS 84 / UTM zone 45N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",87],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32645"]]

    GeoTransform :

    619410.0 30.0 0.0 3089790.0 0.0 30.0

    array(0)

• Data variables: (3)
  • vv
    (acq_date, y, x)
    float32
    dask.array<chunksize=(1, 396, 290), meta=np.ndarray>

    Array Chunk Bytes 45.12 MiB 448.59 kiB Shape (103, 396, 290) (1, 396, 290) Count 1190 Graph Layers 103 Chunks Type float32 numpy.ndarray

  • vh
    (acq_date, y, x)
    float32
    dask.array<chunksize=(1, 396, 290), meta=np.ndarray>

    Array Chunk Bytes 45.12 MiB 448.59 kiB Shape (103, 396, 290) (1, 396, 290) Count 1190 Graph Layers 103 Chunks Type float32 numpy.ndarray

  • ls
    (acq_date, y, x)
    float32
    dask.array<chunksize=(1, 396, 290), meta=np.ndarray>

    Array Chunk Bytes 45.12 MiB 448.59 kiB Shape (103, 396, 290) (1, 396, 290) Count 1190 Graph Layers 103 Chunks Type float32 numpy.ndarray

• Attributes: (0)

Read in vector data

manual outlines of proglacial lakes

lakes = gpd.read_file('/home/emmamarshall/Desktop/siparcs/s1_book/data/proglacial_lake_outlines/proglacial_lake1.shp')
lakes_prj = lakes.to_crs('EPSG:32645')
lakes_prj

	id	geometry
0	1	POLYGON ((623555.903 3099600.816, 623594.809 3...
1	2	POLYGON ((622405.254 3097521.050, 622268.598 3...

Glacier outlines from Randolph Glacier Inventory

#da_bbox = get_bbox_single(asf_clip, 100)

#rgi = gpd.read_file('/home/emmamarshall/Desktop/data/rgi/south_asia_east_15/15_rgi60_SouthAsiaEast.shp')
#rgi.head(3)
#rgi_prj = rgi.to_crs('epsg:32645')

#rgi_sub = gpd.sjoin(rgi_prj, da_bbox, how='inner')

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

power_to_db(asf_clip.vv.mean(dim=['x','y']).plot(ax=axs[0], cmap=plt.cm.Greys_r)
power_to_db(asf_clip.vv.sel(time='2022-04-30').plot(ax=axs[1], cmap=plt.cm.Greys_r)


#rgi_sub.plot(edgecolor='r', facecolor='none', ax=axs[0])
#rgi_sub.plot(edgecolor='r', facecolor='none', ax=axs[1])

lakes_prj.plot(ax=axs[0], facecolor='none', edgecolor='blue')
lakes_prj.plot(ax=axs[1], facecolor='none', edgecolor='blue')

axs[0].set_title('ASF RTC image')
axs[1].set_title('Planetary Computer RTC image')
fig.suptitle('RTC backscatter images of Bhutan glaciers, 30 Apr 2022', fontsize=14)