Explaining the Entire Process of Satellite Data Acquisition: Mastering openEO, STAC, xarray, and rioxarray [Sentinel-2 Data Access Edition]

Introduction

In introductory books on satellite data analysis, the question of "how to acquire satellite data" often becomes the first hurdle. By grasping the general outline of what the actual program is doing before reading the book's explanations, your understanding when reading the book will be significantly faster.

In this article, I will carefully explain the processing performed by the notebook in Chapter 3, Section 1, "Acquiring Satellite Image Data," along with the reasons why it is done that way.

Technologies covered in this article

$$
\begin{array}{|l|l|}
\hline
\textbf{Library / Service} & \textbf{Role} \\
\hline
\text{CDSE (Copernicus Data Space Ecosystem)} & \text{Platform for Sentinel satellite data provided by ESA} \\
\hline
\text{openEO} & \text{Open API standard for satellite data processing} \\
\hline
\text{xarray} & \text{Manipulation of multi-dimensional arrays (NetCDF support)} \\
\hline
\text{rioxarray} & \text{Extension adding GIS functionality to xarray} \\
\hline
\text{STAC} & \text{Standard specification for satellite data catalogs} \\
\hline
\text{pystac-client} & \text{Python client for STAC catalogs} \\
\hline
\text{rasterio} & \text{Reading and writing raster data (GeoTIFF, etc.)} \\
\hline
\text{GDAL} & \text{Standard library for geospatial data conversion} \\
\hline
\end{array}
$$

Overall program structure (flow)

This notebook acquires Sentinel-2 data through roughly two routes.

衛星データ取得
├── ルート1: CDSE + openEO（本書メインルート）
│   ├── 1. CDSEバックエンドに接続
│   ├── 2. OIDC認証
│   ├── 3. DataCube構築（検索条件設定）
│   ├── 4. NetCDF形式でダウンロード
│   ├── 5. xarrayで読み込み・操作
│   ├── 6. SCL（雲マスク用）の可視化
│   ├── 7. RGBトゥルーカラー画像の表示
│   └── 8. rioxarrayでGeoTIFF出力
│
└── ルート2: STAC + pystac-client（参考ルート）
    ├── 1. STACカタログへの接続・検索
    ├── 2. GeoDataFrameで結果整理・雲量ソート
    ├── 3. サムネイル画像で目視確認
    ├── 4. COGファイルのダウンロード
    ├── 5. rasterioでGeoTIFF作成・マスク処理
    └── 6. GDALで8bit変換・表示用画像作成

Route 1: Data acquisition via CDSE + openEO

Why use CDSE?

In October 2023, the old platform for Sentinel satellite data, "Copernicus Open Access Hub," was shut down. Consequently, data acquisition using the previously widely used sentinelsat library is no longer possible.

On the new platform, CDSE (Copernicus Data Space Ecosystem), the following data is available for free (as of April 2024).

• Sentinel-1 / 2 / 3 / 5P / 6 (SAR, optical, atmospheric, etc.)

• Landsat-5 / 7 / 8

• Copernicus DEM (Digital Elevation Model)

• MODIS (Terra / Aqua)

• Many others

User registration should be performed according to the official documentation. Please note that email address verification is required.

What is openEO?

openEO (open Earth Observation) is an international standard for searching, processing, and downloading satellite data using a unified API. A key feature is that you can access multiple platforms with the same code without being conscious of the backend (server-side processing environment).

In the program, we connect to the CDSE backend using openeo.connect("openeo.dataspace.copernicus.eu").

Connection and Authentication Mechanism

Connecting to the Backend

After connecting, calling connection.describe_collection("SENTINEL2_L2A") returns the metadata for that collection (dataset). This allows you to verify that the connection was successful. The returned information includes details such as the geographic coverage, available bands, and the temporal range.

OIDC Authentication

Authentication is required to actually acquire data. In the program, we execute connection.authenticate_oidc().

OIDC (OpenID Connect) is an authentication protocol based on OAuth2. When executed, a URL is output; accessing this in a browser and logging in/authorizing with your CDSE account will issue a token. Once authenticated, the token is cached locally, so it will be used automatically from the next time onwards.

The Concept of DataCube and load_collection

The core concept of openEO is the DataCube.

At the point when connection.load_collection() is called, the data is not yet downloaded. Instead, only an instruction set (DataCube) specifying 'which data is desired' is constructed. The actual data transfer occurs when the download command is issued (lazy evaluation).

The search conditions set in the program are as follows.

Why use Level-2A?

Sentinel-2 data includes Level-1C (top-of-atmosphere reflectance without atmospheric correction) and Level-2A (surface reflectance with atmospheric correction). For many applications such as vegetation analysis and land cover classification, Level-2A is standardly used as it removes atmospheric effects.

Meaning of the Bands

Downloading as NetCDF and the Japanese Path Issue

When you download() a DataCube, a NetCDF (.nc) format file is generated.

NetCDF (Network Common Data Form) is a standard format for storing multidimensional array data such as meteorological, oceanographic, and satellite data. It is also the format that xarray excels at handling the most.

There is one important trick in the program.

# Windowsで日本語パスへの保存を避けるため、ホームディレクトリに保存
NC_PATH = os.path.join(os.path.expanduser("~"), "s2-amami.nc")

This is a workaround for the bug where the netCDF4 C library cannot handle Japanese paths on Windows. We use os.path.expanduser("~") to retrieve the path to the home directory, which consists only of English characters (e.g., C:\Users\username).

Furthermore, we include resample_spatial(resolution=10, projection="EPSG:32652") during the download. There are two reasons for this.

1. Standardizing resolution: Sentinel-2 bands have different resolutions of 10m, 20m, and 60m. Standardizing to 10m ensures that pixels between bands align.

2. Explicit CRS specification: This is a workaround for a server-side bug where CRS information for some tiles is corrupted. EPSG:32652 is the UTM coordinate system (Zone 52N), which is a projected coordinate system suitable for the Amami Islands near Japan.

Data manipulation with xarray

When you load the downloaded NetCDF using xarray.open_dataset(), you get a Dataset object.

The structure of this Dataset is as follows.

Dimensions: (t: N, y: M, x: L)

Coordinates:
  * t    (t)    datetime64[ns]   ← 日時（複数シーン分）
  * x    (x)    float64          ← X座標（投影座標）
  * y    (y)    float64          ← Y座標（投影座標）

Data variables:
    crs   → 座標参照系の情報（配列ではない）
    B04   (t, y, x) float32   ← 赤バンド
    B03   (t, y, x) float32   ← 緑バンド
    B02   (t, y, x) float32   ← 青バンド
    SCL   (t, y, x) float32   ← Scene Classification Layer

Slicing along the time dimension

By specifying the date and time as a string, such as ds.sel(t="2024-03-11"), you can extract data for only that day. This is because xarray supports label-based selection in a way similar to Pandas index operations.

Conversion to RGB data

To combine multiple bands into a single DataArray, use to_array(dim="bands"). This results in a 4D array with dimensions (bands, t, y, x). You can then plot it directly using plot.imshow().

The vmin=0, vmax=2000 settings for display correspond to the scale of DN (Digital Number) values for Sentinel-2 Level-2A. Raw values range from approximately 0 to 10,000, and using the 0 to 2,000 range results in a rendering that is slightly dark to properly exposed.

What is the SCL (Scene Classification Layer)?

The SCL is an auxiliary layer that classifies what each pixel represents. In the program, we display them chronologically using ds["SCL"].plot.imshow(col='t').

Cloud-related values are 3, 8, 9, and 10. In the analysis performed in later chapters, we will use this SCL to mask cloud pixels.

Saving to GeoTIFF with rioxarray

rioxarray is an extension library that adds functionality to xarray for handling spatial reference information (CRS and transformation matrices). You can open a NetCDF with rxr.open_rasterio() and export it to GeoTIFF format using .rio.to_raster().

Note that we specify decode_times=False in open_rasterio(). This is to handle the difference in how the time dimension of NetCDF is decoded between xarray and rioxarray; after reading with decode_times=False, we overwrite ds["t"] with the correctly decoded time information from xarray.

Route 2: Data acquisition via STAC

What is STAC?

STAC (SpatioTemporal Asset Catalog) is a JSON-based specification that standardizes catalog information for satellite data. By publishing STAC-compliant catalogs, data providers allow users to search for and acquire data using a unified API.

The program uses the STAC catalog on AWS at https://earth-search.aws.element84.com/v1. This provides Sentinel-2 data in Cloud-Optimized GeoTIFF (COG) format via the cloud.

Searching and Organizing Results

After connecting to the catalog with pystac_client.Client.open(), you can perform a search using client.search() with the following conditions:

• collections: sentinel-2-l2a (Sentinel-2 Level-2A)

• bbox (Bounding Box): The geographic area to search

• datetime: The date range to search

• query: {"eo:cloud_cover": {"lt": 30}} → Cloud cover less than 30%

Search results are converted into a GeoDataFrame using GeoDataFrame.from_features(). This allows metadata such as cloud cover (eo:cloud_cover) to be treated as columns, making it easy to sort by the least cloudy results using sort_values().

AOI (Area of Interest) Coordinate Transformation

The program may output negative longitude values (such as -220°); this is due to processing that normalizes values exceeding 360°.

# -220は140（東経140°）に変換される
if AREA[i][0] >= 0:
    AREA[i][0] = AREA[i][0] % 360
else:
    AREA[i][0] = -(abs(AREA[i][0]) % 360) + 360

This processing is included because longitude representation methods can vary depending on the data source.

What is a Cloud-Optimized GeoTIFF (COG)?

Much of the data available via STAC is in COG (Cloud-Optimized GeoTIFF) format. Unlike standard GeoTIFFs, COGs are optimized to allow retrieving only the necessary parts via HTTP requests without downloading the entire file. This enables fast acquisition of small areas of interest even from large files.

The URL obtained in the program via selected_item[0][band].href is the URL of the COG file on cloud storage (AWS S3). rasterio supports COGs, allowing you to stream and read files directly from the cloud by passing the URL.

Creating GeoTIFFs and Masking with rasterio

Creating a GeoTIFF

B02 (Blue), B03 (Green), and B04 (Red) are downloaded separately and combined into a single GeoTIFF. When creating a new file with rasterio, the following metadata must be explicitly specified:

• driver: 'Gtiff' (GeoTIFF format)

• width / height: Image size (same as B04)

• count: Number of bands (3 in this case)

• crs: Projected coordinate system (same as B04, epsg:32654)

• transform: Affine transformation matrix (mapping between pixel coordinates and geographic coordinates)

• dtype: Data type (uint16: integers from 0 to 65535)

The reason epsg:32654 (UTM Zone 54N) is used for the CRS is that this search area is near the Japanese mainland (around the Kanto region).

Masking with AOI

Use rasterio.mask.mask() to crop only the polygon area specified. At this time, if the image CRS and the polygon CRS are not matched, an error will occur. In the program, we apply .to_crs(crs='epsg:32654') to the GeoDataFrame to align the coordinate systems.

8-bit conversion using GDAL

Sentinel-2 DN values are uint16 (16-bit), but 8-bit (0-255) is required for standard display. We convert them using the GDAL Translate() function.

-scale 0 255 0 25

This scale parameter means 'convert the original 0-255 to an output of 0-25'. Since actual Sentinel-2 DN values are around 0-10000, this is a very bright setting. You need to adjust the scale according to your purpose (e.g., -scale 0 3000 0 255, etc.).

Comparison of each route

Summary of this article

Next steps

After mastering this notebook, you will learn the following in the subsequent chapters.

• ch3.2: Handling coordinate systems (CRS conversion, projection transformation)

• ch3.3: Detailed usage of GDAL

• ch4: Band arithmetic (NDVI, etc.), analysis of vegetation, roads, farmland, and coastlines

• ch5: Satellite data analysis using machine learning (linear regression, SVM)

• ch6: Unsupervised classification (clustering)

By thoroughly understanding the "data entry point" of satellite data acquisition, it becomes easier to grasp the overall picture before moving on to the next analysis steps.

This article is intended to provide a technical explanation of the sample program (https://github.com/tamanome/satelliteBook) which is published under the Apache License 2.0.