Sentinel-1 data in Python¶

OpenGeoHub Summer School 2023

  • Lorena Abad
  • 2023-08-31

Querying S1-SLC data¶

Sentinel-1 data comes at different levels and provides different products. For applications such as measuring deformation due to tectonic or volcanic activity, quantifying ground subsidence or to generate digital elevation models (DEM), interferometric SAR (InSAR) techniques can be used.

To apply such workflows with Sentinel data, we can use Sentinel-1 Level 1 Single Look Complex products.

SM mode is designed to support ERS (European Remote Sensing) and Envisat missions; IW mode is the default mode over land; EW mode is designed for maritime, ice, and polar zone observation services where wide coverage and short revisit times are demanded; and WV mode is the default mode over the open ocean.

So far, very few cloud computing capabilities are available to compute such complex workflows, therefore, there is still a need to download data. Depending on the application, we will need to download data with certain characteristics.

InSAR principles

Figure from: Xiong S, Muller J-P, Li G. The Application of ALOS/PALSAR InSAR to Measure Subsurface Penetration Depths in Deserts. Remote Sensing. 2017; 9(6):638. https://doi.org/10.3390/rs9060638

Note on terminology.

For DEM generation, for example, we would require a pair of Sentinel-1 scenes acquired closely in time and that have a perpendicular baseline between 150 and 300 m. Usually, computing the perpendicular baseline between two images requires the download of the image pairs.

To avoid downloading several unnecessary Sentinel-1 scenes, we can make use of the Alaska Satellite Facility (ASF) geographic and baseline tools to query the data we need via their API.

Libraries needed for this exercise are imported below:

In [1]:
import asf_search as asf
import geopandas as gpd
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

Define extent¶

We will define an aoi and a start and end date for our queries.

In [2]:
aoi = gpd.read_file("../../data/poznan.geojson")
aoi.explore()
Out[2]:
Make this Notebook Trusted to load map: File -> Trust Notebook
In [7]:
footprint = aoi.to_wkt()
date_start = "2022/05/01"
date_end = "2022/10/01"
footprint.geometry[0]
Out[7]:
'POLYGON ((16.757384 52.535648, 16.757384 52.293644, 17.163617 52.293644, 17.163617 52.535648, 16.757384 52.535648))'

Geographical search¶

Now we can use the asf_search Python module to perform our geographical search. We specify here the platform and the processing level (SLC) that we are looking for, and we limit the results for this exercise to 10 scenes.

In [15]:
products = asf.geo_search(platform=[asf.PLATFORM.SENTINEL1],
                          intersectsWith=footprint.geometry[0],
                          processingLevel=[asf.PRODUCT_TYPE.SLC],
                          beamSwath='IW',
                          start=date_start,
                          end=date_end,
                          maxResults=10)

We can then add the results of the query to a pandas dataframe for easier inspection:

In [14]:
products
Out[14]:
ASFSearchResults([])
In [10]:
products_df = pd.DataFrame([p.properties for p in products])
products_df
Out[10]:
beamModeType browse bytes centerLat centerLon faradayRotation fileID flightDirection groupID granuleType ... processingDate processingLevel sceneName sensor startTime stopTime url pgeVersion fileName frameNumber
0 IW None 4657751122 52.2492 17.9661 None S1A_IW_SLC__1SDV_20220929T163603_20220929T1636... ASCENDING S1A_IWDV_0167_0173_045222_175 SENTINEL_1A_FRAME ... 2022-09-29T16:36:03.000Z SLC S1A_IW_SLC__1SDV_20220929T163603_20220929T1636... C-SAR 2022-09-29T16:36:03.000Z 2022-09-29T16:36:30.000Z https://datapool.asf.alaska.edu/SLC/SA/S1A_IW_... 003.52 S1A_IW_SLC__1SDV_20220929T163603_20220929T1636... 167
1 IW None 4580763215 52.6218 18.4544 None S1A_IW_SLC__1SDV_20220926T050057_20220926T0501... DESCENDING S1A_IWDV_0416_0423_045171_124 SENTINEL_1A_FRAME ... 2022-09-26T05:00:57.000Z SLC S1A_IW_SLC__1SDV_20220926T050057_20220926T0501... C-SAR 2022-09-26T05:00:57.000Z 2022-09-26T05:01:24.000Z https://datapool.asf.alaska.edu/SLC/SA/S1A_IW_... 003.52 S1A_IW_SLC__1SDV_20220926T050057_20220926T0501... 417
2 IW None 4609579705 53.0143 15.6840 None S1A_IW_SLC__1SDV_20220922T164428_20220922T1644... ASCENDING S1A_IWDV_0169_0176_045120_073 SENTINEL_1A_FRAME ... 2022-09-22T16:44:28.000Z SLC S1A_IW_SLC__1SDV_20220922T164428_20220922T1644... C-SAR 2022-09-22T16:44:28.000Z 2022-09-22T16:44:56.000Z https://datapool.asf.alaska.edu/SLC/SA/S1A_IW_... 003.52 S1A_IW_SLC__1SDV_20220922T164428_20220922T1644... 170
3 IW None 4646654586 51.5302 16.1188 None S1A_IW_SLC__1SDV_20220922T164404_20220922T1644... ASCENDING S1A_IWDV_0165_0170_045120_073 SENTINEL_1A_FRAME ... 2022-09-22T16:44:04.000Z SLC S1A_IW_SLC__1SDV_20220922T164404_20220922T1644... C-SAR 2022-09-22T16:44:04.000Z 2022-09-22T16:44:31.000Z https://datapool.asf.alaska.edu/SLC/SA/S1A_IW_... 003.52 S1A_IW_SLC__1SDV_20220922T164404_20220922T1644... 165
4 IW None 4686861613 53.0102 16.5167 None S1A_IW_SLC__1SDV_20220919T050904_20220919T0509... DESCENDING S1A_IWDV_0415_0421_045069_022 SENTINEL_1A_FRAME ... 2022-09-19T05:09:04.000Z SLC S1A_IW_SLC__1SDV_20220919T050904_20220919T0509... C-SAR 2022-09-19T05:09:04.000Z 2022-09-19T05:09:31.000Z https://datapool.asf.alaska.edu/SLC/SA/S1A_IW_... 003.52 S1A_IW_SLC__1SDV_20220919T050904_20220919T0509... 415
5 IW None 4652270299 52.2492 17.9662 None S1A_IW_SLC__1SDV_20220917T163602_20220917T1636... ASCENDING S1A_IWDV_0167_0173_045047_175 SENTINEL_1A_FRAME ... 2022-09-17T16:36:02.000Z SLC S1A_IW_SLC__1SDV_20220917T163602_20220917T1636... C-SAR 2022-09-17T16:36:02.000Z 2022-09-17T16:36:30.000Z https://datapool.asf.alaska.edu/SLC/SA/S1A_IW_... 003.52 S1A_IW_SLC__1SDV_20220917T163602_20220917T1636... 167
6 IW None 4734573721 52.6220 18.4523 None S1A_IW_SLC__1SDV_20220914T050057_20220914T0501... DESCENDING S1A_IWDV_0416_0423_044996_124 SENTINEL_1A_FRAME ... 2022-09-14T05:00:57.000Z SLC S1A_IW_SLC__1SDV_20220914T050057_20220914T0501... C-SAR 2022-09-14T05:00:57.000Z 2022-09-14T05:01:24.000Z https://datapool.asf.alaska.edu/SLC/SA/S1A_IW_... 003.52 S1A_IW_SLC__1SDV_20220914T050057_20220914T0501... 417
7 IW None 4707081403 53.0145 15.6814 None S1A_IW_SLC__1SDV_20220910T164429_20220910T1644... ASCENDING S1A_IWDV_0169_0176_044945_073 SENTINEL_1A_FRAME ... 2022-09-10T16:44:29.000Z SLC S1A_IW_SLC__1SDV_20220910T164429_20220910T1644... C-SAR 2022-09-10T16:44:29.000Z 2022-09-10T16:44:56.000Z https://datapool.asf.alaska.edu/SLC/SA/S1A_IW_... 003.52 S1A_IW_SLC__1SDV_20220910T164429_20220910T1644... 170
8 IW None 4721237800 51.5303 16.1161 None S1A_IW_SLC__1SDV_20220910T164404_20220910T1644... ASCENDING S1A_IWDV_0164_0171_044945_073 SENTINEL_1A_FRAME ... 2022-09-10T16:44:04.000Z SLC S1A_IW_SLC__1SDV_20220910T164404_20220910T1644... C-SAR 2022-09-10T16:44:04.000Z 2022-09-10T16:44:32.000Z https://datapool.asf.alaska.edu/SLC/SA/S1A_IW_... 003.52 S1A_IW_SLC__1SDV_20220910T164404_20220910T1644... 165
9 IW None 4575317676 53.0107 16.5132 None S1A_IW_SLC__1SDV_20220907T050904_20220907T0509... DESCENDING S1A_IWDV_0415_0421_044894_022 SENTINEL_1A_FRAME ... 2022-09-07T05:09:04.000Z SLC S1A_IW_SLC__1SDV_20220907T050904_20220907T0509... C-SAR 2022-09-07T05:09:04.000Z 2022-09-07T05:09:31.000Z https://datapool.asf.alaska.edu/SLC/SA/S1A_IW_... 003.52 S1A_IW_SLC__1SDV_20220907T050904_20220907T0509... 415

10 rows × 28 columns

Ascending vs. Descending:

Schema pass - Ascending pass - Descending pass"/>

Baseline search¶

Now that we have scenes that intersect with our defined extent, we can do a baseline search that will allow us to fetch all the S1 scenes that pair with the first S1 result from our geographical query. The baseline search returns a set of products with precomputed perpendicular baselines, so that we can focus our download on the data that we need.

In [16]:
stack = products[0].stack()
In [17]:
print(f'{len(stack)} products found in stack')

446 products found in stack

We can take a look at the data again as a pandas data frame and we will see that the last two columns correspond to the temporal and perpendicular baseline.

In [19]:
stack_df = pd.DataFrame([p.properties for p in stack])
stack_df.sort_values('perpendicularBaseline')
Out[19]:
beamModeType browse bytes centerLat centerLon faradayRotation fileID flightDirection groupID granuleType ... sceneName sensor startTime stopTime url pgeVersion fileName frameNumber temporalBaseline perpendicularBaseline
416 IW None 4563266039 52.2489 17.9638 None S1A_IW_SLC__1SDV_20220905T163603_20220905T1636... ASCENDING S1A_IWDV_0167_0173_044872_175 SENTINEL_1A_FRAME ... S1A_IW_SLC__1SDV_20220905T163603_20220905T1636... C-SAR 2022-09-05T16:36:03.000Z 2022-09-05T16:36:30.000Z https://datapool.asf.alaska.edu/SLC/SA/S1A_IW_... 003.52 S1A_IW_SLC__1SDV_20220905T163603_20220905T1636... 167 -24 -184.0
111 IW None 4745376346 51.4741 18.1578 None S1B_IW_SLC__1SDV_20170504T163429_20170504T1634... ASCENDING S1B_IWDV_0165_0170_005451_175 SENTINEL_1B_FRAME ... S1B_IW_SLC__1SDV_20170504T163429_20170504T1634... C-SAR 2017-05-04T16:34:29.000Z 2017-05-04T16:34:56.000Z https://datapool.asf.alaska.edu/SLC/SB/S1B_IW_... 002.82 S1B_IW_SLC__1SDV_20170504T163429_20170504T1634... 165 -1974 -117.0
110 IW None 5629299887 53.0986 17.6801 None S1B_IW_SLC__1SDV_20170504T163454_20170504T1635... ASCENDING S1B_IWDV_0170_0176_005451_175 SENTINEL_1B_FRAME ... S1B_IW_SLC__1SDV_20170504T163454_20170504T1635... C-SAR 2017-05-04T16:34:54.000Z 2017-05-04T16:35:26.000Z https://datapool.asf.alaska.edu/SLC/SB/S1B_IW_... 002.82 S1B_IW_SLC__1SDV_20170504T163454_20170504T1635... 170 -1974 -114.0
321 IW None 4696301380 52.2496 17.9639 None S1A_IW_SLC__1SDV_20201009T163551_20201009T1636... ASCENDING S1A_IWDV_0167_0173_034722_175 SENTINEL_1A_FRAME ... S1A_IW_SLC__1SDV_20201009T163551_20201009T1636... C-SAR 2020-10-09T16:35:51.000Z 2020-10-09T16:36:18.000Z https://datapool.asf.alaska.edu/SLC/SA/S1A_IW_... 003.31 S1A_IW_SLC__1SDV_20201009T163551_20201009T1636... 167 -720 -99.0
99 IW None 4698539592 51.4743 18.1576 None S1B_IW_SLC__1SDV_20170317T163428_20170317T1634... ASCENDING S1B_IWDV_0165_0170_004751_175 SENTINEL_1B_FRAME ... S1B_IW_SLC__1SDV_20170317T163428_20170317T1634... C-SAR 2017-03-17T16:34:28.000Z 2017-03-17T16:34:54.000Z https://datapool.asf.alaska.edu/SLC/SB/S1B_IW_... 002.72 S1B_IW_SLC__1SDV_20170317T163428_20170317T1634... 165 -2022 -91.0
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
219 IW None 4664583363 52.2503 17.9663 None S1A_IW_SLC__1SDV_20190205T163535_20190205T1636... ASCENDING S1A_IWDV_0167_0173_025797_175 SENTINEL_1A_FRAME ... S1A_IW_SLC__1SDV_20190205T163535_20190205T1636... C-SAR 2019-02-05T16:35:35.000Z 2019-02-05T16:36:02.000Z https://datapool.asf.alaska.edu/SLC/SA/S1A_IW_... 002.91 S1A_IW_SLC__1SDV_20190205T163535_20190205T1636... 167 -1332 173.0
429 IW None 4536281700 52.2498 17.9679 None S1A_IW_SLC__1SDV_20230208T163559_20230208T1636... ASCENDING S1A_IWDV_0167_0173_047147_175 SENTINEL_1A_FRAME ... S1A_IW_SLC__1SDV_20230208T163559_20230208T1636... C-SAR 2023-02-08T16:35:59.000Z 2023-02-08T16:36:26.000Z https://datapool.asf.alaska.edu/SLC/SA/S1A_IW_... 003.52 S1A_IW_SLC__1SDV_20230208T163559_20230208T1636... 167 132 178.0
72 IW None 6584072949 53.2640 17.6334 None S1B_IW_SLC__1SDV_20161129T163455_20161129T1635... ASCENDING S1B_IWDV_0170_0177_003176_175 SENTINEL_1B_FRAME ... S1B_IW_SLC__1SDV_20161129T163455_20161129T1635... C-SAR 2016-11-29T16:34:55.000Z 2016-11-29T16:35:33.000Z https://datapool.asf.alaska.edu/SLC/SB/S1B_IW_... 002.72 S1B_IW_SLC__1SDV_20161129T163455_20161129T1635... 170 -2130 185.0
425 IW None 4732915777 52.2497 17.9671 None S1A_IW_SLC__1SDV_20221222T163601_20221222T1636... ASCENDING S1A_IWDV_0167_0173_046447_175 SENTINEL_1A_FRAME ... S1A_IW_SLC__1SDV_20221222T163601_20221222T1636... C-SAR 2022-12-22T16:36:01.000Z 2022-12-22T16:36:28.000Z https://datapool.asf.alaska.edu/SLC/SA/S1A_IW_... 003.52 S1A_IW_SLC__1SDV_20221222T163601_20221222T1636... 167 84 191.0
9 IW None 4887202715 52.5500 17.8553 None S1A_IW_SLC__1SDV_20150310T163515_20150310T1635... ASCENDING S1A_IWDV_0168_0173_004972_175 SENTINEL_1A_FRAME ... S1A_IW_SLC__1SDV_20150310T163515_20150310T1635... C-SAR 2015-03-10T16:35:15.000Z 2015-03-10T16:35:43.000Z https://datapool.asf.alaska.edu/SLC/SA/S1A_IW_... 002.36 S1A_IW_SLC__1SDV_20150310T163515_20150310T1635... 168 -2760 NaN

446 rows × 30 columns

To have an idea of how spread our data is, we can plot the temporal and the perpendicular baselines against each other.

In [20]:
stack_df.plot.scatter(x="temporalBaseline", y="perpendicularBaseline")
Out[20]:
<Axes: xlabel='temporalBaseline', ylabel='perpendicularBaseline'>
No description has been provided for this image

Ideally, we will filter those values where temporalBaseline <= 30 and 150 <= perpendicularBaseline <= 300 for instance to get image pairs suitable for DEM generation. So we can filter our data frame for those values. We look for absolute values since the order of the images is not relevant.

In [21]:
stack_df[(abs(stack_df['temporalBaseline']) <= 30) &
         (abs(stack_df['perpendicularBaseline']) >= 150) &
         (abs(stack_df['perpendicularBaseline']) <= 300)]
Out[21]:
beamModeType browse bytes centerLat centerLon faradayRotation fileID flightDirection groupID granuleType ... sceneName sensor startTime stopTime url pgeVersion fileName frameNumber temporalBaseline perpendicularBaseline
416 IW None 4563266039 52.2489 17.9638 None S1A_IW_SLC__1SDV_20220905T163603_20220905T1636... ASCENDING S1A_IWDV_0167_0173_044872_175 SENTINEL_1A_FRAME ... S1A_IW_SLC__1SDV_20220905T163603_20220905T1636... C-SAR 2022-09-05T16:36:03.000Z 2022-09-05T16:36:30.000Z https://datapool.asf.alaska.edu/SLC/SA/S1A_IW_... 003.52 S1A_IW_SLC__1SDV_20220905T163603_20220905T1636... 167 -24 -184.0

1 rows × 30 columns

We only get one image fitting the characteristics we require. Let's look at its properties:

In [22]:
stack[416].properties
Out[22]:
{'beamModeType': 'IW',
 'browse': None,
 'bytes': 4563266039,
 'centerLat': 52.2489,
 'centerLon': 17.9638,
 'faradayRotation': None,
 'fileID': 'S1A_IW_SLC__1SDV_20220905T163603_20220905T163630_044872_055BFF_5F88-SLC',
 'flightDirection': 'ASCENDING',
 'groupID': 'S1A_IWDV_0167_0173_044872_175',
 'granuleType': 'SENTINEL_1A_FRAME',
 'insarStackId': None,
 'md5sum': '72e8baa2456d336cb9b71e2c5f7e93be',
 'offNadirAngle': None,
 'orbit': 44872,
 'pathNumber': 175,
 'platform': 'Sentinel-1A',
 'pointingAngle': None,
 'polarization': 'VV+VH',
 'processingDate': '2022-09-05T16:36:03.000Z',
 'processingLevel': 'SLC',
 'sceneName': 'S1A_IW_SLC__1SDV_20220905T163603_20220905T163630_044872_055BFF_5F88',
 'sensor': 'C-SAR',
 'startTime': '2022-09-05T16:36:03.000Z',
 'stopTime': '2022-09-05T16:36:30.000Z',
 'url': 'https://datapool.asf.alaska.edu/SLC/SA/S1A_IW_SLC__1SDV_20220905T163603_20220905T163630_044872_055BFF_5F88.zip',
 'pgeVersion': '003.52',
 'fileName': 'S1A_IW_SLC__1SDV_20220905T163603_20220905T163630_044872_055BFF_5F88.zip',
 'frameNumber': 167,
 'temporalBaseline': -24,
 'perpendicularBaseline': -184}

Let's also remember this is paired with the original product we calcualted the baselines for.

In [23]:
products[0].properties
Out[23]:
{'beamModeType': 'IW',
 'browse': None,
 'bytes': 4657751122,
 'centerLat': 52.2492,
 'centerLon': 17.9661,
 'faradayRotation': None,
 'fileID': 'S1A_IW_SLC__1SDV_20220929T163603_20220929T163630_045222_0567B7_46CF-SLC',
 'flightDirection': 'ASCENDING',
 'groupID': 'S1A_IWDV_0167_0173_045222_175',
 'granuleType': 'SENTINEL_1A_FRAME',
 'insarStackId': None,
 'md5sum': 'a37a41c9155ed7a93b4055066710e224',
 'offNadirAngle': None,
 'orbit': 45222,
 'pathNumber': 175,
 'platform': 'Sentinel-1A',
 'pointingAngle': None,
 'polarization': 'VV+VH',
 'processingDate': '2022-09-29T16:36:03.000Z',
 'processingLevel': 'SLC',
 'sceneName': 'S1A_IW_SLC__1SDV_20220929T163603_20220929T163630_045222_0567B7_46CF',
 'sensor': 'C-SAR',
 'startTime': '2022-09-29T16:36:03.000Z',
 'stopTime': '2022-09-29T16:36:30.000Z',
 'url': 'https://datapool.asf.alaska.edu/SLC/SA/S1A_IW_SLC__1SDV_20220929T163603_20220929T163630_045222_0567B7_46CF.zip',
 'pgeVersion': '003.52',
 'fileName': 'S1A_IW_SLC__1SDV_20220929T163603_20220929T163630_045222_0567B7_46CF.zip',
 'frameNumber': 167}

Downloading the data¶

Finally, with the ASF API we can download our data to further analyse it with, e.g. SNAP. To do so we can make use of the url property.

In [ ]:
urls = [
    products[0].properties['url'],
    stack[416].properties['url']
]

Once that is set we can use the download_urls() function as speccified below to get our data in a desired directory. To download the data we will need EarthData credentials. This notebook from the ASF describes the authentication process.

asf.download_urls(urls=urls, path='data/s1', session=user_pass_session, processes=5)

Exploring S1-RTC data¶

In [24]:
import matplotlib.pyplot as plt
import numpy as np
import os
import pandas as pd
import rioxarray as rio
import xarray as xr

Now let's take a look at a bit more processed data that we can directly work with. Still in Level-1 you will see the Ground Range Detected (GRD) product in the figure above. This is S1 data that has been further processed (it has been detected, multi-looked and projected to ground range). The SLC products we queried before preserve phase information and are processed at the natural pixel spacing whereas GRD products contain the detected amplitude and are multi-looked to reduce the impact of speckle.

An extra processing step is to perform Radiometric Terrain Correction, and some data providers like Microsoft Planetary Computer make this dataset available worldwide. Feel free to explore the Planetary Computer access options to work on larger datasets if you are interested.

In the spirit to avoid the need for you to get credentials for this particular workshop, we will use a Sentinel-1 RTC dataset for the Contiguous United States (CONUS) which is freely accessible.

We will access this data using the Amazon Web Services (AWS) CLI directly (with the awscli package). Let's explore the available data:

In [31]:
!aws s3 ls s3://sentinel-s1-rtc-indigo/tiles/RTC/1/IW/10/S/DH/2016/S1A_20160729_10SDH_ASC/ --no-sign-request

                           PRE S1A_IW_GRDH_1SSV_20160729T020719_20160729T020744_012357_013406_0D26/
2020-09-25 13:22:52   32864309 Gamma0_VV.tif
2020-09-25 13:23:01    1351178 local_incident_angle.tif

To download a scene we can directly request the data as:

In [ ]:
# Run only if you want to have 110MB of random data on your laptop!
!aws s3 cp s3://sentinel-s1-rtc-indigo/tiles/RTC/1/IW/12/R/UV/2021/S1B_20210121_12RUV_DSC/Gamma0_VV.tif S1B_20210121_12RUV_DSC/Gamma0_VV.tif --no-sign-request

The available bands have the prefix Gamma0. This is the result of the RTC algorithm. Read more about the backscatter types here.

We will also see that the data has a suffix, either VV or VH, this is the polarization. That refers to the way data is collected.

Polarization HH-VV

Read more about it here, here and here.

Let's start exploring the data. For this we will use rioxarray. We will set an environment key to establish no sign request for AWS. And we will also be leveraging the tight integration between xarray and dask to lazily read in data via the chunks parameter.

We will get both the vv and vh data.

In [32]:
os.environ['AWS_NO_SIGN_REQUEST'] = 'YES'
In [33]:
url_vv = 's3://sentinel-s1-rtc-indigo/tiles/RTC/1/IW/10/U/CU/2019/S1B_20190109_10UCU_ASC/Gamma0_VV.tif'
url_vh = 's3://sentinel-s1-rtc-indigo/tiles/RTC/1/IW/10/U/CU/2019/S1B_20190109_10UCU_ASC/Gamma0_VH.tif'
s1_vv = rio.open_rasterio(url_vv, chunks=True)
s1_vh = rio.open_rasterio(url_vh, chunks=True)
In [34]:
s1_vv
Out[34]:
<xarray.DataArray (band: 1, y: 5490, x: 5490)>
dask.array<open_rasterio-19cbe1dffda706ec0d650050fad8aee1<this-array>, shape=(1, 5490, 5490), dtype=float32, chunksize=(1, 5490, 5490), chunktype=numpy.ndarray>
Coordinates:
  * band         (band) int64 1
  * x            (x) float64 3e+05 3e+05 3e+05 ... 4.098e+05 4.098e+05 4.098e+05
  * y            (y) float64 5.4e+06 5.4e+06 5.4e+06 ... 5.29e+06 5.29e+06
    spatial_ref  int64 0
Attributes: (12/18)
    ABSOLUTE_ORBIT_NUMBER:  14411
    AREA_OR_POINT:          Area
    DATE:                   2019-01-09
    MISSION_ID:             S1B
    NUMBER_SCENES:          2
    ORBIT_DIRECTION:        ascending
    ...                     ...
    TILE_ID:                10UCU
    VALID_PIXEL_PERCENT:    87.946
    _FillValue:             0.0
    scale_factor:           1.0
    add_offset:             0.0
    long_name:              Gamma0_VV
xarray.DataArray
  • band: 1
  • y: 5490
  • x: 5490
  • dask.array<chunksize=(1, 5490, 5490), meta=np.ndarray>
    Array Chunk
    Bytes 114.98 MiB 114.98 MiB
    Shape (1, 5490, 5490) (1, 5490, 5490)
    Dask graph 1 chunks in 2 graph layers
    Data type float32 numpy.ndarray
    		5490 5490 1
    • band
      (band)
      int64
      1
      array([1])
    • x
      (x)
      float64
      3e+05 3e+05 ... 4.098e+05 4.098e+05
      array([300010., 300030., 300050., ..., 409750., 409770., 409790.])
    • y
      (y)
      float64
      5.4e+06 5.4e+06 ... 5.29e+06
      array([5399990., 5399970., 5399950., ..., 5290250., 5290230., 5290210.])
    • spatial_ref
      ()
      int64
      0
      crs_wkt :
      PROJCS["WGS 84 / UTM zone 10N",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",-123],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","32610"]]
      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 10N
      grid_mapping_name :
      transverse_mercator
      latitude_of_projection_origin :
      0.0
      longitude_of_central_meridian :
      -123.0
      false_easting :
      500000.0
      false_northing :
      0.0
      scale_factor_at_central_meridian :
      0.9996
      spatial_ref :
      PROJCS["WGS 84 / UTM zone 10N",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",-123],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","32610"]]
      GeoTransform :
      300000.0 20.0 0.0 5400000.0 0.0 -20.0
      array(0)
    • band
      PandasIndex
      PandasIndex(Index([1], dtype='int64', name='band'))
    • x
      PandasIndex
      PandasIndex(Index([300010.0, 300030.0, 300050.0, 300070.0, 300090.0, 300110.0, 300130.0,
       300150.0, 300170.0, 300190.0,
       ...
       409610.0, 409630.0, 409650.0, 409670.0, 409690.0, 409710.0, 409730.0,
       409750.0, 409770.0, 409790.0],
      dtype='float64', name='x', length=5490))
    • y
      PandasIndex
      PandasIndex(Index([5399990.0, 5399970.0, 5399950.0, 5399930.0, 5399910.0, 5399890.0,
       5399870.0, 5399850.0, 5399830.0, 5399810.0,
       ...
       5290390.0, 5290370.0, 5290350.0, 5290330.0, 5290310.0, 5290290.0,
       5290270.0, 5290250.0, 5290230.0, 5290210.0],
      dtype='float64', name='y', length=5490))
  • ABSOLUTE_ORBIT_NUMBER :
    14411
    AREA_OR_POINT :
    Area
    DATE :
    2019-01-09
    MISSION_ID :
    S1B
    NUMBER_SCENES :
    2
    ORBIT_DIRECTION :
    ascending
    OVR_RESAMPLING_ALG :
    AVERAGE
    SCENES :
    S1B_IW_GRDH_1SDV_20190109T020931_20190109T020956_014411_01AD36_7CA8,S1B_IW_GRDH_1SDV_20190109T020956_20190109T021021_014411_01AD36_9989
    SCENE_1_METADATA :
    {"title": "S1B_IW_GRDH_1SDV_20190109T020931_20190109T020956_014411_01AD36_7CA8", "mission_id": "S1B", "sensor_operational_mode": "IW", "product_type": "RTC", "resolution_class": "H", "processing_level": "1", "polarization_mode": "DV", "start_time": "2019-01-09T02:09:31", "end_time": "2019-01-09T02:09:56", "absolute_orbit_number": "14411", "mission_data_take_id": "01AD36", "product_unique_identifier": "7CA8", "footprint": "{\"type\": \"Polygon\", \"coordinates\": [[[-126.535973, 47.750252], [-123.125053, 48.149876], [-122.782219, 46.652382], [-126.098656, 46.253361], [-126.535973, 47.750252]]]}", "orbit_direction": "ASC", "ingested_at": "2019-01-09T07:16:25.694000", "size": "914219674"}
    SCENE_1_PRODUCT_INFO :
    {"id": "S1B_IW_GRDH_1SDV_20190109T020931_20190109T020956_014411_01AD36_7CA8", "path": "scenes/RTC/1/2019/1/9/IW/DV/S1B_IW_GRDH_1SDV_20190109T020931_20190109T020956_014411_01AD36_7CA8/", "missionId": "S1B", "productType": "RTC", "mode": "IW", "polarization": "DV", "startTime": "2019-01-09T02:09:31", "stopTime": "2019-01-09T02:09:56", "absoluteOrbitNumber": 14411, "missionDataTakeId": "01AD36", "productUniqueIdentifier": "7CA8", "sciHubIngestion": "2019-01-09T07:16:25.694000Z", "s3Ingestion": "2020-09-23T18:01:22.339964Z", "sciHubId": "073d0c7d-47c7-44dd-be5f-bea7479c3cfa", "footprint": {"type": "Polygon", "coordinates": [[[-126.535973, 47.750252], [-123.125053, 48.149876], [-122.782219, 46.652382], [-126.098656, 46.253361], [-126.535973, 47.750252]]]}}
    SCENE_2_METADATA :
    {"title": "S1B_IW_GRDH_1SDV_20190109T020956_20190109T021021_014411_01AD36_9989", "mission_id": "S1B", "sensor_operational_mode": "IW", "product_type": "RTC", "resolution_class": "H", "processing_level": "1", "polarization_mode": "DV", "start_time": "2019-01-09T02:09:56", "end_time": "2019-01-09T02:10:21", "absolute_orbit_number": "14411", "mission_data_take_id": "01AD36", "product_unique_identifier": "9989", "footprint": "{\"type\": \"Polygon\", \"coordinates\": [[[-126.984451, 49.24667], [-123.469872, 49.647362], [-123.123749, 48.150101], [-126.535995, 47.75034], [-126.984451, 49.24667]]]}", "orbit_direction": "ASC", "ingested_at": "2019-01-09T07:40:51.331000", "size": "981747202"}
    SCENE_2_PRODUCT_INFO :
    {"id": "S1B_IW_GRDH_1SDV_20190109T020956_20190109T021021_014411_01AD36_9989", "path": "scenes/RTC/1/2019/1/9/IW/DV/S1B_IW_GRDH_1SDV_20190109T020956_20190109T021021_014411_01AD36_9989/", "missionId": "S1B", "productType": "RTC", "mode": "IW", "polarization": "DV", "startTime": "2019-01-09T02:09:56", "stopTime": "2019-01-09T02:10:21", "absoluteOrbitNumber": 14411, "missionDataTakeId": "01AD36", "productUniqueIdentifier": "9989", "sciHubIngestion": "2019-01-09T07:40:51.331000Z", "s3Ingestion": "2020-09-23T18:01:22.340412Z", "sciHubId": "709311d1-d033-487d-875a-786ff0570ebc", "footprint": {"type": "Polygon", "coordinates": [[[-126.984451, 49.24667], [-123.469872, 49.647362], [-123.123749, 48.150101], [-126.535995, 47.75034], [-126.984451, 49.24667]]]}}
    TILE_ID :
    10UCU
    VALID_PIXEL_PERCENT :
    87.946
    _FillValue :
    0.0
    scale_factor :
    1.0
    add_offset :
    0.0
    long_name :
    Gamma0_VV

Let's visualize a slice of the data:

In [35]:
s1_vv_ss = s1_vv.isel(x=slice(5000, 5500), y=slice(1000, 1500)).compute()
In [36]:
s1_vv_ss.plot(cmap=plt.cm.Greys_r)
Out[36]:
<matplotlib.collections.QuadMesh at 0x7f19112ea950>
No description has been provided for this image

To better visualize the data, we can apply a power to dB scale. This transformation applies a logarithmic scale to the data for easier visualisation, but it is not recommended to use this for any computations, since the data gets distorted.

In [37]:
def power_to_db(input_arr):
    return (10*np.log10(np.abs(input_arr)))
In [38]:
power_to_db(s1_vv_ss).plot(cmap=plt.cm.Greys_r)
Out[38]:
<matplotlib.collections.QuadMesh at 0x7f191136b490>
No description has been provided for this image

Exercise:¶

  • Try to combine the s1_vv and s1_vh objects together, compute a new band with the result of VH/VV and use these three layers to generate a false color RGB composite.
In [39]:
s1_rt = s1_vh/s1_vv
s1 = xr.concat([s1_vh, s1_vv, s1_rt], pd.Index(['vh', 'vv', 'vhvv'], name = 'band'))

s1
Out[39]:
<xarray.DataArray (band: 3, y: 5490, x: 5490)>
dask.array<concatenate, shape=(3, 5490, 5490), dtype=float32, chunksize=(1, 5490, 5490), chunktype=numpy.ndarray>
Coordinates:
  * band         (band) object 'vh' 'vv' 'vhvv'
  * x            (x) float64 3e+05 3e+05 3e+05 ... 4.098e+05 4.098e+05 4.098e+05
  * y            (y) float64 5.4e+06 5.4e+06 5.4e+06 ... 5.29e+06 5.29e+06
    spatial_ref  int64 0
Attributes: (12/18)
    ABSOLUTE_ORBIT_NUMBER:  14411
    AREA_OR_POINT:          Area
    DATE:                   2019-01-09
    MISSION_ID:             S1B
    NUMBER_SCENES:          2
    ORBIT_DIRECTION:        ascending
    ...                     ...
    TILE_ID:                10UCU
    VALID_PIXEL_PERCENT:    87.946
    _FillValue:             0.0
    scale_factor:           1.0
    add_offset:             0.0
    long_name:              Gamma0_VH
xarray.DataArray
  • band: 3
  • y: 5490
  • x: 5490
  • dask.array<chunksize=(1, 5490, 5490), meta=np.ndarray>
    Array Chunk
    Bytes 344.93 MiB 114.98 MiB
    Shape (3, 5490, 5490) (1, 5490, 5490)
    Dask graph 3 chunks in 6 graph layers
    Data type float32 numpy.ndarray
    		5490 5490 3
    • band
      (band)
      object
      'vh' 'vv' 'vhvv'
      array(['vh', 'vv', 'vhvv'], dtype=object)
    • x
      (x)
      float64
      3e+05 3e+05 ... 4.098e+05 4.098e+05
      array([300010., 300030., 300050., ..., 409750., 409770., 409790.])
    • y
      (y)
      float64
      5.4e+06 5.4e+06 ... 5.29e+06
      array([5399990., 5399970., 5399950., ..., 5290250., 5290230., 5290210.])
    • spatial_ref
      ()
      int64
      0
      crs_wkt :
      PROJCS["WGS 84 / UTM zone 10N",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",-123],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","32610"]]
      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 10N
      grid_mapping_name :
      transverse_mercator
      latitude_of_projection_origin :
      0.0
      longitude_of_central_meridian :
      -123.0
      false_easting :
      500000.0
      false_northing :
      0.0
      scale_factor_at_central_meridian :
      0.9996
      spatial_ref :
      PROJCS["WGS 84 / UTM zone 10N",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",-123],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","32610"]]
      GeoTransform :
      300000.0 20.0 0.0 5400000.0 0.0 -20.0
      array(0)
    • band
      PandasIndex
      PandasIndex(Index(['vh', 'vv', 'vhvv'], dtype='object', name='band'))
    • x
      PandasIndex
      PandasIndex(Index([300010.0, 300030.0, 300050.0, 300070.0, 300090.0, 300110.0, 300130.0,
       300150.0, 300170.0, 300190.0,
       ...
       409610.0, 409630.0, 409650.0, 409670.0, 409690.0, 409710.0, 409730.0,
       409750.0, 409770.0, 409790.0],
      dtype='float64', name='x', length=5490))
    • y
      PandasIndex
      PandasIndex(Index([5399990.0, 5399970.0, 5399950.0, 5399930.0, 5399910.0, 5399890.0,
       5399870.0, 5399850.0, 5399830.0, 5399810.0,
       ...
       5290390.0, 5290370.0, 5290350.0, 5290330.0, 5290310.0, 5290290.0,
       5290270.0, 5290250.0, 5290230.0, 5290210.0],
      dtype='float64', name='y', length=5490))
  • ABSOLUTE_ORBIT_NUMBER :
    14411
    AREA_OR_POINT :
    Area
    DATE :
    2019-01-09
    MISSION_ID :
    S1B
    NUMBER_SCENES :
    2
    ORBIT_DIRECTION :
    ascending
    OVR_RESAMPLING_ALG :
    AVERAGE
    SCENES :
    S1B_IW_GRDH_1SDV_20190109T020931_20190109T020956_014411_01AD36_7CA8,S1B_IW_GRDH_1SDV_20190109T020956_20190109T021021_014411_01AD36_9989
    SCENE_1_METADATA :
    {"title": "S1B_IW_GRDH_1SDV_20190109T020931_20190109T020956_014411_01AD36_7CA8", "mission_id": "S1B", "sensor_operational_mode": "IW", "product_type": "RTC", "resolution_class": "H", "processing_level": "1", "polarization_mode": "DV", "start_time": "2019-01-09T02:09:31", "end_time": "2019-01-09T02:09:56", "absolute_orbit_number": "14411", "mission_data_take_id": "01AD36", "product_unique_identifier": "7CA8", "footprint": "{\"type\": \"Polygon\", \"coordinates\": [[[-126.535973, 47.750252], [-123.125053, 48.149876], [-122.782219, 46.652382], [-126.098656, 46.253361], [-126.535973, 47.750252]]]}", "orbit_direction": "ASC", "ingested_at": "2019-01-09T07:16:25.694000", "size": "914219674"}
    SCENE_1_PRODUCT_INFO :
    {"id": "S1B_IW_GRDH_1SDV_20190109T020931_20190109T020956_014411_01AD36_7CA8", "path": "scenes/RTC/1/2019/1/9/IW/DV/S1B_IW_GRDH_1SDV_20190109T020931_20190109T020956_014411_01AD36_7CA8/", "missionId": "S1B", "productType": "RTC", "mode": "IW", "polarization": "DV", "startTime": "2019-01-09T02:09:31", "stopTime": "2019-01-09T02:09:56", "absoluteOrbitNumber": 14411, "missionDataTakeId": "01AD36", "productUniqueIdentifier": "7CA8", "sciHubIngestion": "2019-01-09T07:16:25.694000Z", "s3Ingestion": "2020-09-23T18:00:51.876503Z", "sciHubId": "073d0c7d-47c7-44dd-be5f-bea7479c3cfa", "footprint": {"type": "Polygon", "coordinates": [[[-126.535973, 47.750252], [-123.125053, 48.149876], [-122.782219, 46.652382], [-126.098656, 46.253361], [-126.535973, 47.750252]]]}}
    SCENE_2_METADATA :
    {"title": "S1B_IW_GRDH_1SDV_20190109T020956_20190109T021021_014411_01AD36_9989", "mission_id": "S1B", "sensor_operational_mode": "IW", "product_type": "RTC", "resolution_class": "H", "processing_level": "1", "polarization_mode": "DV", "start_time": "2019-01-09T02:09:56", "end_time": "2019-01-09T02:10:21", "absolute_orbit_number": "14411", "mission_data_take_id": "01AD36", "product_unique_identifier": "9989", "footprint": "{\"type\": \"Polygon\", \"coordinates\": [[[-126.984451, 49.24667], [-123.469872, 49.647362], [-123.123749, 48.150101], [-126.535995, 47.75034], [-126.984451, 49.24667]]]}", "orbit_direction": "ASC", "ingested_at": "2019-01-09T07:40:51.331000", "size": "981747202"}
    SCENE_2_PRODUCT_INFO :
    {"id": "S1B_IW_GRDH_1SDV_20190109T020956_20190109T021021_014411_01AD36_9989", "path": "scenes/RTC/1/2019/1/9/IW/DV/S1B_IW_GRDH_1SDV_20190109T020956_20190109T021021_014411_01AD36_9989/", "missionId": "S1B", "productType": "RTC", "mode": "IW", "polarization": "DV", "startTime": "2019-01-09T02:09:56", "stopTime": "2019-01-09T02:10:21", "absoluteOrbitNumber": 14411, "missionDataTakeId": "01AD36", "productUniqueIdentifier": "9989", "sciHubIngestion": "2019-01-09T07:40:51.331000Z", "s3Ingestion": "2020-09-23T18:00:51.876961Z", "sciHubId": "709311d1-d033-487d-875a-786ff0570ebc", "footprint": {"type": "Polygon", "coordinates": [[[-126.984451, 49.24667], [-123.469872, 49.647362], [-123.123749, 48.150101], [-126.535995, 47.75034], [-126.984451, 49.24667]]]}}
    TILE_ID :
    10UCU
    VALID_PIXEL_PERCENT :
    87.946
    _FillValue :
    0.0
    scale_factor :
    1.0
    add_offset :
    0.0
    long_name :
    Gamma0_VH
In [40]:
s1_ss = s1.isel(x=slice(5000, 5500), y=slice(1000, 1500)).compute()

/opt/conda/lib/python3.11/site-packages/dask/core.py:121: RuntimeWarning: invalid value encountered in divide
  return func(*(_execute_task(a, cache) for a in args))
In [41]:
power_to_db(s1_ss.sel(band=["vv", "vh", "vhvv"])).plot.imshow(rgb="band", robust=True)
Out[41]:
<matplotlib.image.AxesImage at 0x7f19106c6950>
No description has been provided for this image

More resources:¶

Working with Sentinel-1 SLC data can also be done with Python. There are a couple of packages available for this (snappy and snapista), but ESA is currently working on a follow up of the snappy package called esa-snappy which will be compatible with the upcoming SNAP-10. Since the developers claim its worth the wait, I would at this point direct you to the webpage where they seem to be documenting basic usage of the tool. So for that feel free to check this site more or less at the end of August.

For more complex InSAR analysis on the cloud, there are some other alternatives. Take a look at this paper by Piter et al. (2021) for an overview and comparison. There are some other attempts still under development, e.g., dinoSAR.

A lot of the examples for this notebook, mainly for the RTC processing were adapted from Emma Marshall's excellent tutorial on Sentinel-1 RTC data workflows with xarray.