{
"cells": [
{
"cell_type": "markdown",
"id": "brutal-tutorial",
"metadata": {},
"source": [
"# Sentinel-1\n",
"\n",
"```{admonition} Learning Objectives\n",
"*A 30 minute guide to Sentinel-1 data for SnowEX*\n",
"- understand key characteristics of Sentinel-1 Synthetic Aperture Radar\n",
"- find, visualize, interpret Sentinel-1 data products\n",
"- use Python raster libraries [rioxarray](https://corteva.github.io/rioxarray) and [hvplot](https://hvplot.holoviz.org)\n",
"```\n",
"\n",
":::{figure-md} sentinel1\n",
"\n",
"\n",
"Artist view of Sentinel-1. image source: `https://www.esa.int/ESA_Multimedia/Images/2014/01/Sentinel-1_radar_vision`)\n",
":::\n",
"\n",
"```{seealso}\n",
"this tutorial is a quick practical guide and *will not cover InSAR processing*, check out [UNAVCO InSAR Short Courses](https://www.unavco.org/education/professional-development/short-courses/short-courses.html) if you're interested in learning custom processing of SAR data.\n",
"```"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "classified-namibia",
"metadata": {
"tags": [
"remove-input"
]
},
"outputs": [],
"source": [
"# Import all Python libraries required for this notebook\n",
"import hvplot.xarray\n",
"import os\n",
"import pandas as pd\n",
"import rioxarray\n",
"import s3fs"
]
},
{
"cell_type": "markdown",
"id": "double-tucson",
"metadata": {},
"source": [
"## Dive right in\n",
"\n",
"Synthetic Aperture Radar (SAR) is an active imaging technique that records microwave reflections off of Earth's surface. **Unlike passive optical imagers that require cloud-free, sunny days, SAR can operate at night and microwaves penetrate clouds** At first glance, a SAR 'image' might look a lot like a black-and-white image of the Earth, but these observations contain more than just color values and can be used to query many physical properties and processes! \n",
"\n",
"But before getting into theory and caveats, let's visualize some data. An easy way to get started with Sentinel-1 SAR over a SnowEx field site is to use the Radiometric Terrain Corrected (RTC) backscatter amplitude data on AWS: https://registry.opendata.aws/sentinel-1-rtc-indigo/ "
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "honest-active",
"metadata": {},
"outputs": [],
"source": [
"# GDAL environment variables to efficiently read remote data\n",
"os.environ['GDAL_DISABLE_READDIR_ON_OPEN']='EMPTY_DIR' \n",
"os.environ['AWS_NO_SIGN_REQUEST']='YES' \n",
"\n",
"# Data is stored in a public S3 Bucket\n",
"url = 's3://sentinel-s1-rtc-indigo/tiles/RTC/1/IW/12/S/YJ/2016/S1B_20161121_12SYJ_ASC/Gamma0_VV.tif'\n",
"\n",
"# These Cloud-Optimized-Geotiff (COG) files have 'overviews', low-resolution copies for quick visualization\n",
"da = rioxarray.open_rasterio(url, overview_level=3).squeeze('band')\n",
"da.hvplot.image(clim=(0,0.4), cmap='gray', \n",
" x='x', y='y', \n",
" aspect='equal', frame_width=400,\n",
" title='S1B_20161121_12SYJ_ASC',\n",
" rasterize=True # send rendered image to browser, rather than full array\n",
" )"
]
},
{
"cell_type": "markdown",
"id": "smooth-disclosure",
"metadata": {},
"source": [
"```{admonition} Interpretation\n",
"The above image is in UTM coordinates, with linear power units. Dark zones (such as Grand Mesa) correspond to low radar amplitude returns, which can result from wet snow. High-amplitude bright zones occur on dry slopes that perpendicular to radar incidence and urban areas like Grand Junction, CO.\n",
"```"
]
},
{
"cell_type": "markdown",
"id": "inclusive-programmer",
"metadata": {},
"source": [
"```{admonition} Exercise\n",
":class: dropdown\n",
"\n",
"Visualizations of SAR data are often better on a different scale. \n",
"Convert the linear power values to amplitude or decibels and replot.\n",
"Need a hint? Check out this [short article](https://storymaps.arcgis.com/stories/73b6af082e1f44bca8a0c5fb6bf09f37)\n",
"```"
]
},
{
"cell_type": "markdown",
"id": "rotary-extent",
"metadata": {},
"source": [
"## Quick facts\n",
"\n",
"Sentinel-1 is a constellation of two C-band satellites operated by the European Space Agency (ESA). *It is the first SAR system with a global monitoring strategy and fully open data policy!* S1A launched 3 April 2014, and S1B launched 25 April 2016. There are many observation modes for SAR, over land the most common mode is 'Interferometric Wideswath' (IW), which has the following characteristics:\n",
"\n",
"| wavelength | resolution | posting | frame size | incidence | orbit repeat |\n",
"| - | - | - | - | - | - | \n",
"| *(cm)* | *rng x azi (m)* | *rng x azi (m)* | *w x h (km)* | *(deg)* | *(days)* |\n",
"| 5.547 | 3x22 | 2.3x14 | 250x170 | 33-43 | 12 |\n",
"\n",
"\n",
"Unlike most optical satellite observations, SAR antennas are pointed to the side, resulting in an \"line-of-sight\" (LOS) incidence angle with respect to the ellipsoid normal vector. A consequence of this viewing geometry is that radar images can have significant distortions known as shadow (for example due to tall mountains) and layover (where large regions perpendicular to incident energy all map to the same pixel). Also note that *resolution* != *pixel posting*. *resolution* is the minimal separate of distinguishable targets on the ground, and *posting* is the raster grid of a Sentinel1 image. \n",
"\n",
"The [global observation plan](https://sentinels.copernicus.eu/web/sentinel/missions/sentinel-1/observation-scenario) for Sentinel-1 has changed over time, but generally over Europe, you have repeat acquisitions every 6 days, and elsewhere every 12 or 24 days.\n",
"\n",
"Satellite radar instruments record the *amplitude* and *phase* of microwaves that bounce off Earth's surface. These values can be stored as a single complex number, so you'll often encounter SAR images that store complex-valued arrays. {term}`InSAR`, or 'Interferometric Synthetic Aperture Radar', is a common processing pipline to generate phase-change maps or 'interferograms' (which can be related to tectonic surface displacements or snow properties) using two different SAR acquisitions.\n",
"\n",
":::{figure-md} sar-schematic\n",
"\n",
"\n",
"InSAR schematic depicting how phase shifts can result from earthquake displacements. [Source](https://www.ga.gov.au/scientific-topics/positioning-navigation/geodesy/geodetic-techniques/interferometric-synthetic-aperture-radar))\n",
":::\n",
"\n",
"```\n",
"\n",
"```{seealso}\n",
"- [ESA Documentation](https://sentinel.esa.int/web/sentinel/missions/sentinel-1)\n",
"- [Alaska Satellite Facility Documentation](https://asf.alaska.edu/data-sets/sar-data-sets/sentinel-1/)\n",
"```"
]
},
{
"cell_type": "markdown",
"id": "described-twenty",
"metadata": {},
"source": [
"## Search and Discovery\n",
"\n",
"### Public Archives\n",
"The most common data products you'll encounter for Sentinel-1 are {term}`GRD` (just amplitude) and {term}`SLC` (amplitude+phase). These are [level-1](https://earthdata.nasa.gov/collaborate/open-data-services-and-software/data-information-policy/data-levels) products that many higher-level scientific products are derived from. \n",
"\n",
"Level-1 data is typically stored in *radar* coordinates. For SLCs, pixels are rectangular with a ratio of about 6:1 for range:azimuth, given the different resolutions in these orthogonal directions. You might notice ground control points (GCPs) stored in GRD and SLC metadata for [approximate geocoding](https://docs.qgis.org/3.16/en/docs/user_manual/working_with_raster/georeferencer.html#available-transformation-algorithms), but you need to use advanced processing pipelines to accurately geocode and generate higher level products such as interferograms {cite:p}`YagueMartinez2016`.\n",
"\n",
"### Higher-level products\n",
"There are also a growing number of cloud-hosted archives of processed Sentinel-1 data including global $\\sigma_0$ radiometric terrain corrected (RTC) in [Google Earth Engine](https://developers.google.com/earth-engine/guides/sentinel1), $\\gamma_0$ RTC data as an [AWS public dataset](https://registry.opendata.aws/sentinel-1-rtc-indigo/). Generating higher-level products often requires using a digital elevation model, so it's important to be aware of the *resolution* and *acquisition time* of that digital elevation model. For example, two common free global digital elevation models for processing include:\n",
"\n",
"| Product| Acquisition Date | Coverage | Resolution (m) | Open Access |\n",
"| - | - | - | - | - | \n",
"| [SRTM](https://lpdaac.usgs.gov/products/srtmgl1v003/) | 2000-02-11 | 60$^\\circ$N to 60$^\\circ$S | 30 | [link](https://lpdaac.usgs.gov/products/srtmgl1v003/) |\n",
"| [Copernicus DEM](https://portal.opentopography.org/datasetMetadata?otCollectionID=OT.032021.4326.1) | 2019-12-01 -> | global | 30 | [link](https://registry.opendata.aws/copernicus-dem/) | \n",
"\n",
"\n",
"### On-demand processing\n",
"SAR processing algorithms are complicated and resource intensive, so when possible it is nice to take advantages of services to generate higher level products. For example ASF has the HyP3 service which can generate RTC and Interferometric products: https://hyp3-docs.asf.alaska.edu"
]
},
{
"cell_type": "markdown",
"id": "pressed-strength",
"metadata": {},
"source": [
"## Amplitude\n",
"\n",
"SAR Backscatter variations can be related to of melting snow {cite:p}`Small2011`. Cross-polarized Sentinel-1 backscatter variation has even been shown to relate to {term}`SWE` {cite:p}`Lievens2019`. Let's use the public RTC data from earlier and interpret values over grand mesa.\n",
"\n",
"In order to work with a multidimensional timeseries stack of imagery we'll construct a [GDAL VRT file](https://gdal.org/drivers/raster/vrt.html) Based on the following organization: `s3://sentinel-s1-rtc-indigo/tiles/RTC/1/[MODE]/[MGRS UTM zone]/[MGRS latitude label]/[MGRS Grid Square ID]/[YEAR]/[SATELLITE]_[DATE]_[TILE ID]_[ORBIT DIRECTION]/[ASSET]`, since the code takes a few minutes to run, we only run it if the vrt doesn't already exist:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "forced-multimedia",
"metadata": {
"tags": [
"hide-input"
]
},
"outputs": [],
"source": [
"zone = 12\n",
"latLabel = 'S'\n",
"square = 'YJ'\n",
"year = '202*' #>=2020\n",
"date = '*' #all acquisitions\n",
"polarization = 'VV'\n",
"s3Path = f's3://sentinel-s1-rtc-indigo/tiles/RTC/1/IW/{zone}/{latLabel}/{square}/{year}/{date}/Gamma0_{polarization}.tif'\n",
"\n",
"# Find imagery according to S3 path pattern\n",
"s3 = s3fs.S3FileSystem(anon=True)\n",
"keys = s3.glob(s3Path[5:]) #strip s3://\n",
"print(f'Located {len(keys)} images matching {s3Path}:')\n",
"\n",
"vrtName = f'stack{zone}{latLabel}{square}.vrt'\n",
"if not os.path.exists(vrtName):\n",
" with open('s3paths.txt', 'w') as f:\n",
" for key in keys:\n",
" f.write(\"/vsis3/%s\\n\" % key)\n",
"\n",
" cmd = f'gdalbuildvrt -overwrite -separate -input_file_list s3paths.txt {vrtName}'\n",
" print(cmd)\n",
" os.system(cmd)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "twelve-property",
"metadata": {},
"outputs": [],
"source": [
"# Load a time series we created a VRT with GDAL to facilitate this step\n",
"da3 = rioxarray.open_rasterio(vrtName, overview_level=3, chunks='auto')\n",
"da3; #omit output"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "fuzzy-johnson",
"metadata": {},
"outputs": [],
"source": [
"# Need to add time coordinates to this data\n",
"datetimes = [pd.to_datetime(x[55:63]) for x in keys]\n",
" \n",
"# add new coordinate to existing dimension \n",
"da = da3.assign_coords(time=('band', datetimes))\n",
"# make 'time' active coordinate instead of integer band\n",
"da = da.swap_dims({'band':'time'})\n",
"# Name the dataset (helpful for hvplot calls later on)\n",
"da.name = 'Gamma0VV'\n",
"da"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "executive-saturday",
"metadata": {},
"outputs": [],
"source": [
"#use a small bounding box over grand mesa (UTM coordinates)\n",
"xmin,xmax,ymin,ymax = [739186, 742748, 4.325443e+06, 4.327356e+06]\n",
"daT = da.sel(x=slice(xmin, xmax), \n",
" y=slice(ymax, ymin))\n",
"\n",
"# NOTE: this can take a while on slow internet connections, we're reading over 100 images!\n",
"all_points = daT.where(daT!=0).hvplot.scatter('time', groupby=[], dynspread=True, datashade=True) \n",
"mean_trend = daT.where(daT!=0, drop=True).mean(dim=['x','y']).hvplot.line(title='North Grand Mesa', color='red')\n",
"(all_points * mean_trend)"
]
},
{
"cell_type": "markdown",
"id": "insured-necklace",
"metadata": {},
"source": [
"```{admonition} Interpretation\n",
"The background backscatter for this area of interest is approximately 0.5. A distinct dip in backscatter is observed during Spring snow melt April through May, with max backscatter in June corresponding to bare ground conditions. Decreasing backscatter amplitude from June onwards is likely due to increasing vegetation. \n",
"```"
]
},
{
"cell_type": "markdown",
"id": "aware-blackjack",
"metadata": {},
"source": [
"## Phase\n",
"\n",
"InSAR phase delays can be due to a number of factors including tectonic motions, atmospheric water vapor changes, and ionospheric conditions. The theory relating InSAR phase delays due to propagation in dry snow is described in a classic study by {cite:p}`Guneriussen2001`.\n",
"\n",
"A first-order approximation of Snow-Water-Equivalent from InSAR phase change from {cite}`Leinss2015` is:\n",
"\n",
"$$\n",
" \\Delta SWE = \\frac{\\Delta \\Phi \\lambda_i}{2 \\pi (1.59 + \\theta_i^{5/2})}\n",
"$$ (phase_swe) \n",
"\n",
"Where $\\Delta\\Phi$ is measured change in phase, $\\lambda_i$ is the radar wavelength and $\\theta_i$ is the incidence angle. This approximation assumes dry, homogeneous snow with a depth of less than 3 meters. Note also that phase delays are also be caused by changes in atmospheric water vapor, ionospheric conditions, and tectonic displacements, so care must be taken to isolate phase changes arising from SWE changes. Isolating these signals is complicated and more studies like SnowEx are necessary to validate satellite-based SWE extractions with in-situ sensors.\n",
"\n",
"The following cell gets you started with plotting phase data generated by ASF's on-demand InSAR processor. It takes about an hour for processing an interferogram, so we've done that ahead of time with the scripts in this repository https://github.com/snowex-hackweek/hyp3SAR, and data outputs here https://github.com/snowex-hackweek/tutorial-data."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "dirty-render",
"metadata": {},
"outputs": [],
"source": [
"%%bash \n",
"\n",
"# Retrieve a copy of data files used in this tutorial from Zenodo.org:\n",
"# Re-running this cell will not re-download things if they already exist\n",
"\n",
"mkdir -p /tmp/tutorial-data\n",
"cd /tmp/tutorial-data\n",
"wget -q -nc -O data.zip https://zenodo.org/record/5504396/files/sar.zip\n",
"unzip -q -n data.zip\n",
"rm data.zip"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "suitable-referral",
"metadata": {},
"outputs": [],
"source": [
"path = '/tmp/tutorial-data/sar/sentinel1/S1AA_20201030T131820_20201111T131820_VVP012_INT80_G_ueF_EBD2/S1AA_20201030T131820_20201111T131820_VVP012_INT80_G_ueF_EBD2_unw_phase.tif'\n",
"da = rioxarray.open_rasterio(path, masked=True).squeeze('band')"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "bound-resort",
"metadata": {},
"outputs": [],
"source": [
"da.hvplot.image(x='x', y='y', aspect='equal', rasterize=True, cmap='plasma', title='2020/10/30_2020/11/11 Unwrapped Phase (radians)')"
]
},
{
"cell_type": "markdown",
"id": "fantastic-bacon",
"metadata": {},
"source": [
"```{admonition} Interpretation\n",
":class: danger\n",
"Single Sentinel-1 interferograms are usually dominated by atmospheric noise, so be cautious about interpretations of surface properties such as surface displacements or SWE changes. In theory, positive changes correspond to a path increases (or propagation delay). Phase changes are expected to be minimal for short durations, so it's common to apply a DC-offset so that the image has a mean of 0, or normalize the image to a ground-control point with assumed zero phase shift.\n",
"```"
]
},
{
"cell_type": "markdown",
"id": "saved-seventh",
"metadata": {},
"source": [
"```{admonition} Exercises\n",
":class: dropdown\n",
"\n",
"- Plot a histogram of phase change\n",
"- Convert the phase change into a map of SWE change\n",
"- Process another interferogram between separate dates for the same area and average the two\n",
"```"
]
},
{
"cell_type": "markdown",
"id": "signal-letter",
"metadata": {},
"source": [
"## Next steps\n",
"\n",
"This tutorial just scratched the surface of Sentinel-1 SAR! Hopefully you're now eager to explore more and utilize this dataset for SnowEx projects. Below are additional references and project ideas:\n",
"\n",
"```{seealso}\n",
"- expanded tutorial on [AWS RTC public dataset](https://github.com/scottyhq/sentinel1-rtc)\n",
"- documentation for open source [ISCE2](https://github.com/isce-framework/isce2-docs) SAR processing software\n",
"- [UAF Microwave Remote Sensing Course](https://radar.community.uaf.edu)\n",
"- [SERVIR SAR Handbook](https://servirglobal.net/Global/Articles/Article/2674/sar-handbook-comprehensive-methodologies-for-forest-monitoring-and-biomass-estimation) \n",
"```\n",
"\n",
"```{admonition} Project Ideas\n",
":class: tip\n",
"- Explore RTC backscatter over different areas of interest and time ranges\n",
"- Apply the C-Snow algorithm for SWE retrieval from cross-polarized backscatter with the AWS RTC public dataset \n",
"- Process interferograms during dry-snow conditions and convert time series of phase changes due to snow accumulation. This would require looking at historical weather and other in-situ sensors\n",
"- Compare Sentinel-1 and UAVSAR amplitude and phase products (different wavelengths and viewing geometries)\n",
"```"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.8"
}
},
"nbformat": 4,
"nbformat_minor": 5
}