🐛

Docker + ESA SNAP Python で Sentinel-1 をオルソ幾何補正する

2024/01/11に公開

Docker

Python

sar

tech

概要

dockerを用いてESA SNAP Pythonの環境を作成し、Sentinel-1 GRDデータをオルソ幾何補正（Geometric Terrain Correction）する。

環境

Ubuntu 20.04.5 LTS
Docker 20.10.18
QGIS 3.10.4

環境構築

docker build

ESA SNAP 9 docker imagesを元に、以下のコマンドを実行して環境を構築する。このビルドは1-2時間かかる。ここで既存のmundialis/esa-snap:ubuntuのイメージを使用しない理由は、既存のイメージではNo valid orbit file found for 01-JAN-2024 20:42:51.625682 Orbit files may be downloaded from https://scihub.copernicus.eu/gnss/odata/v1/の問題に遭遇し、Access to orbit files from 1 November 2023等の最新のアップデートに対応できていないためである。そのため、SNAPのversionが8であるinstallation_ESA_snappyのイメージも利用できない。

$ git clone https://github.com/mundialis/esa-snap.git -b ubuntu
$ cd esa-snap
$ docker build -t esa-python .

docker run

以下のコマンドを実行する。また今後必要なライブラリをインストールする。

$ docker run -it esa-python:latest 
root@703a13e9f2a3:/# pip install pygeoif pyshp

Sentinel-1 GRD のダウンロード

Copernicus BrowserからSentinel-1 GRDをダウンロードする。今回は以下の画像の東京周辺のデータとして、S1A_IW_GRDH_1SDV_20240101T204356_20240101T204421_051918_0645D3_9D6B.SAFE.zipをダウンロードした。

オルソ幾何補正（Geometric Terrain Correction）の適用

クリップ領域の作成

【QGIS】ポリゴンを作成するを参考にして、ディズニーランド周辺のshapefileを作成し、disney.shpというファイルで保存する。

docker cp

S1A_IW_GRDH_1SDV_20240101T204356_20240101T204421_051918_0645D3_9D6B.SAFE.zip, disney.shp, clip.py, terrain_correction.pyをコンテナ内にコピーする。ここで、clip.py, terrain_correction.pyはFlood mapping using Snappyを参考に以下に記述している。

$ docker cp S1A_IW_GRDH_1SDV_20240101T204356_20240101T204421_051918_0645D3_9D6B.SAFE.zip 703a13e9f2a3:/
$ docker cp disney.shp 703a13e9f2a3:/
$ docker cp clip.py 703a13e9f2a3:/
$ docker cp terrain_correction.py 703a13e9f2a3:/

clip.py

import pygeoif
import shapefile
from snappy import GPF, HashMap, ProductIO, WKTReader, jpy


## Helper functions
def showProductInformation(product):
    width = product.getSceneRasterWidth()
    print("Width: {} px".format(width))
    height = product.getSceneRasterHeight()
    print("Height: {} px".format(height))
    name = product.getName()
    print("Name: {}".format(name))
    band_names = product.getBandNames()
    print("Band names: {}".format(", ".join(band_names)))


def shpToWKT(shp_path):
    r = shapefile.Reader(shp_path)
    g = []
    for s in r.shapes():
        g.append(pygeoif.geometry.as_shape(s))
    m = pygeoif.MultiPoint(g)
    return str(m.wkt).replace("MULTIPOINT", "POLYGON(") + ")"


## Preprocessing functions
def applyOrbit(product):
    parameters = HashMap()
    parameters.put("orbitType", "Sentinel Precise (Auto Download)")
    parameters.put("polyDegree", "3")
    parameters.put("continueOnFail", False)
    return GPF.createProduct("Apply-Orbit-File", parameters, product)


def subset(product, shpPath):
    parameters = HashMap()
    wkt = shpToWKT(shpPath)
    SubsetOp = jpy.get_type("org.esa.snap.core.gpf.common.SubsetOp")
    geometry = WKTReader().read(wkt)
    parameters.put("copyMetadata", True)
    parameters.put("geoRegion", geometry)
    return GPF.createProduct("Subset", parameters, product)


if __name__ == "__main__":
    ## GPF Initialization
    GPF.getDefaultInstance().getOperatorSpiRegistry().loadOperatorSpis()

    file_path = (
        "S1A_IW_GRDH_1SDV_20240101T204356_20240101T204421_051918_0645D3_9D6B.SAFE.zip"
    )
    product = ProductIO.readProduct(file_path)
    showProductInformation(product)

    product_orbitfile = applyOrbit(product)

    product_subset = subset(product_orbitfile, "disney.shp")
    showProductInformation(product_subset)
    ProductIO.writeProduct(product_subset, "product_subset", "GeoTIFF")

terrain_correction.py

from snappy import GPF, HashMap, ProductIO, WKTReader, jpy


## Helper functions
def showProductInformation(product):
    width = product.getSceneRasterWidth()
    print("Width: {} px".format(width))
    height = product.getSceneRasterHeight()
    print("Height: {} px".format(height))
    name = product.getName()
    print("Name: {}".format(name))
    band_names = product.getBandNames()
    print("Band names: {}".format(", ".join(band_names)))


## Preprocessing functions
def calibration(product):
    parameters = HashMap()
    parameters.put("outputSigmaBand", True)
    parameters.put("sourceBands", "Intensity_VV")
    parameters.put("selectedPolarisations", "VV")
    parameters.put("outputImageScaleInDb", False)
    return GPF.createProduct("Calibration", parameters, product)


def speckleFilter(product):
    parameters = HashMap()
    filterSizeY = "5"
    filterSizeX = "5"
    parameters.put("sourceBands", "Sigma0_VV")
    parameters.put("filter", "Lee")
    parameters.put("filterSizeX", filterSizeX)
    parameters.put("filterSizeY", filterSizeY)
    parameters.put("dampingFactor", "2")
    parameters.put("estimateENL", "true")
    parameters.put("enl", "1.0")
    parameters.put("numLooksStr", "1")
    parameters.put("targetWindowSizeStr", "3x3")
    parameters.put("sigmaStr", "0.9")
    parameters.put("anSize", "50")
    return GPF.createProduct("Speckle-Filter", parameters, product)


def terrainCorrection(product):
    parameters = HashMap()
    parameters.put("demName", "SRTM 3Sec")
    parameters.put("pixelSpacingInMeter", 10.0)
    parameters.put("sourceBands", "Sigma0_VV")
    return GPF.createProduct("Terrain-Correction", parameters, product)


if __name__ == "__main__":
    ## GPF Initialization
    GPF.getDefaultInstance().getOperatorSpiRegistry().loadOperatorSpis()

    file_path = "product_subset.tif"
    product = ProductIO.readProduct(file_path)
    showProductInformation(product)

    product_preprocessed = terrainCorrection(speckleFilter(calibration(product)))
    showProductInformation(product_preprocessed)
    ProductIO.writeProduct(product_preprocessed, "product_preprocessed", "GeoTIFF")

クリップ処理

以下のコマンドを実行する。このプログラムはSentinel-1 GRDをdisney.shpでクリップ処理し、その結果をproduct_subset.tifに保存する。

root@703a13e9f2a3:/# python3 clip.py

オルソ幾何補正

以下のコマンドを実行する。このプログラムはproduct_subset.tifにキャリブレーション、スペックルフィルタリング、オルソ幾何補正を適用し、product_preprocessed.tifに保存する。

root@703a13e9f2a3:/# python3 terrain_correction.py

データ確認

product_preprocessed.tifをQGISで可視化すると以下のようになり、クリップ処理、オルソ幾何補正が適用されていることが確認できる。