Skip to main content
  • »
  • EODATA »
  • Integration with Earth Observation Tools (EODATA)

Integration with Earth Observation Tools (EODATA)

Table of Contents

  1. Introduction

  2. Prerequisites

  3. ESA SNAP

  4. GDAL (cmd)

  5. Python Rasterio

Introduction

This section presents practical scenarios for integrating data from the EODATA repository with popular Earth Observation (EO) data processing tools. The examples included cover various levels of complexity — from simple file operations (GDAL), through batch processing in the ESA SNAP environment, to integration with Python using Rasterio.

The goal is to show how EODATA data can be incorporated into the EO (Earth Observation) analysis process depending on needs:

  • quick inspection and conversion,

  • automation of processing chains,

  • integration with the Python analytics ecosystem.

Three groups of tools are used for this purpose: * ESA SNAP — a rich set of algorithms for satellite processing, especially for Sentinel missions. * GDAL — a universal set of command-line tools for raster and vector data operations. * Rasterio — a high-level Python library, ideal for integration with data analysis.

Prerequisites

  • Access to EODATA data (e.g., Sentinel-1/2/3).

  • Installed tools:

    • Docker (for SNAP),

    • GDAL (≥ 3.x),

    • Python 3.x + rasterio.

ESA SNAP

SNAP (Sentinel Application Platform) is a tool dedicated to analyzing data from Sentinel missions. It stands out with:

  • Graph Processing Framework (GPF) — defining processing workflows in XML format,

  • ability to run in batch mode (CLI, Docker),

  • integration with Python through the snappy module.

Example usage in Docker

SNAP can be run in a container, which eliminates configuration issues:

Running SNAP in Docker

docker run -it --rm --name snap esa/snap:latest

Running processing graphs in ESA SNAP

docker run -it --rm mundialis/esa-snap:latest /usr/local/snap/bin/gpt <GraphFile.xml> [options] [<source-file-1> <source-file-2> ...]

Integration with Python (snappy)

import snappy
from snappy import ProductIO, GPF

product = ProductIO.readProduct('/eodata/Sentinel-2/MSI/L2A/2025/09/26/S2C_MSIL2A_20250926T100041_N0511_R122_T34UCD_20250926T152116.SAFE')
parameters = snappy.HashMap()
parameters.put('targetBands', 'B4,B8')
result = GPF.createProduct('Subset', parameters, product)
ProductIO.writeProduct(result, 'subset.dim', 'BEAM-DIMAP')

Automating XML schemas (workflow)

gpt GraphProcessing.xml -Pinput="/eodata/Sentinel-2/MSI/L2A/2025/09/26/S2C_MSIL2A_20250926T100041_N0511_R122_T34UCD_20250926T152116.SAFE" -Poutput="output.dim"

(where GraphProcessing.xml is a file defining the operation workflow, e.g., atmospheric correction, resampling)

GDAL (cmd)

GDAL is a set of command-line tools that allow for quick operations on raster data. It is efficient, easy to automate, and doesn’t require additional dependencies.

Example applications:

Check version:

gdalinfo --version

Display metadata:

gdalinfo /eodata/Sentinel-2/MSI/L2A/2025/09/26/S2C_MSIL2A_20250926T100041_N0511_R122_T34UCD_20250926T152116.SAFE/GRANULE/L2A_T34UCD_A005526_20250926T100041/IMG_DATA/R10m/T34UCD_20250926T100041_B02_10m.jp2

Convert to GeoTIFF:

gdal_translate /eodata/Sentinel-2/MSI/L2A/2025/09/26/S2C_MSIL2A_20250926T100041_N0511_R122_T34UCD_20250926T152116.SAFE/GRANULE/L2A_T34UCD_A005526_20250926T100041/IMG_DATA/R10m/T34UCD_20250926T100041_B02_10m.jp2 output.tif

Reproject to WGS84:

gdalwarp -t_srs EPSG:4326 /eodata/Sentinel-2/MSI/L2A/2025/09/26/S2C_MSIL2A_20250926T100041_N0511_R122_T34UCD_20250926T152116.SAFE/GRANULE/L2A_T34UCD_A005526_20250926T100041/IMG_DATA/R10m/T34UCD_20250926T100041_B02_10m.jp2 T34UCD_20250926T100041_B02_10m_wgs84.tif

Resample to 10m:

gdalwarp -tr 10 10 /eodata/Sentinel-2/MSI/L2A/2025/09/26/S2C_MSIL2A_20250926T100041_N0511_R122_T34UCD_20250926T152116.SAFE/GRANULE/L2A_T34UCD_A005526_20250926T100041/IMG_DATA/R20m/T34UCD_20250926T100041_B12_20m.jp2 T34UCD_20250926T100041_B12_10m_resampled.tif

Clip AOI (based on shapefile):

gdalwarp -cutline aoi.shp -crop_to_cutline /eodata/Sentinel-2/MSI/L2A/2025/09/26/S2C_MSIL2A_20250926T100041_N0511_R122_T34UCD_20250926T152116.SAFE/GRANULE/L2A_T34UCD_A005526_20250926T100041/IMG_DATA/R20m/T34UCD_20250926T100041_B12_20m.jp2 T34UCD_20250926T100041_B12_10m_resampled.jp2 clipped.tif

Calculate NDVI (e.g., Sentinel-2, B8 = NIR, B4 = RED):

gdal_calc.py \
  -A /eodata/Sentinel-2/MSI/L2A/2025/09/26/S2C_MSIL2A_20250926T100041_N0511_R122_T34UCD_20250926T152116.SAFE/GRANULE/L2A_T34UCD_A005526_20250926T100041/IMG_DATA/R10m/T34UCD_20250926T100041_B08_10m.jp2 \
  -B /eodata/Sentinel-2/MSI/L2A/2025/09/26/S2C_MSIL2A_20250926T100041_N0511_R122_T34UCD_20250926T152116.SAFE/GRANULE/L2A_T34UCD_A005526_20250926T100041/IMG_DATA/R10m/T34UCD_20250926T100041_B04_10m.jp2 \
  --outfile=/SSDCUBE/ndvi.tif \
  --calc="(A.astype(float)-B.astype(float))/(A.astype(float)+B.astype(float))" \
  --type=Float32 \
  --NoDataValue=-9999

Compression and optimization (COG):

gdal_translate /eodata/Sentinel-2/MSI/L2A/2025/09/26/S2C_MSIL2A_20250926T100041_N0511_R122_T34UCD_20250926T152116.SAFE/GRANULE/L2A_T34UCD_A005526_20250926T100041/IMG_DATA/R10m/T34UCD_20250926T100041_B04_10m.jp2 T34UCD_20250926T100041_B04_10m_output_cog.tif -of COG

Python Rasterio

Rasterio is a Python library based on GDAL, but offering a convenient high-level interface. It easily integrates with the scientific ecosystem (NumPy, Pandas, GeoPandas, scikit-learn).

Example application:

import rasterio

with rasterio.open("S2A_MSIL1C_20210901_B04.jp2") as src:
    print(src.meta)
    print(src.bounds)

# Convert to GeoTIFF

with rasterio.open("input.jp2") as src:
    profile = src.profile
    with rasterio.open("output.tif", "w", **profile) as dst:
        dst.write(src.read())

# Clip AOI

import geopandas as gpd
from rasterio.mask import mask

shp = gpd.read_file("aoi.shp")
with rasterio.open("input.tif") as src:
    out_image, out_transform = mask(src, shp.geometry, crop=True)

# Calculate NDVI

with rasterio.open("B8.tif") as nir, rasterio.open("B4.tif") as red:
    nir_arr = nir.read(1).astype('float32')
    red_arr = red.read(1).astype('float32')
    ndvi = (nir_arr - red_arr) / (nir_arr + red_arr)