技術資料

衛星画像汀線抽出プログラムCoastSatについて
大塚

淳一

で安価に汀線位置を得ることが可能な衛星画像を用い１．はじめにた汀線抽出技術の開発が進められている4)、5)、6)。Vos
砂浜は背後に控える生活・産業活動の場を高潮、高ら7)、8)が開発したPython tool kit CoastSatはオープン

波から守る防護機能に加え、海水の浄化機能、生態系

衛星画像のLandsatとSentinel-2から高い時間分解能での保全機能、レクリエーションの場を提供する親水機

長期間の汀線（Landsatの場合、過去約40年分、16日能など、多くの優れた機能を有している 。我が国で間隔）を半自動的に抽出することができる。このようは、河川からの土砂供給量の減少や海岸での構造物設な高い時間分解能の長期汀線データを得ることによっ

置に伴う土砂収支の不均衡などの影響により、過去て、高波や津波などの短期的要因のみならず、海面水数十年で砂浜の侵食が進み、その規模は年間約160 ha

位や波候の変化などの長期的要因による汀線変化を定（1978年～1993年の侵食速度）にも及んでいる 。近年、量的に評価することが可能となり、短期・長期両面か侵食対策技術が蓄積され、侵食対策事業により砂浜のら管理手法を検討することができる。また、Landsat保全・回復が図られている一方、対策の多くが侵食にやSentinel-2は今後も運用が続くと予想されることかよる影響が深刻になってから後追い的に行われている

ら、CoastSatは汀線モニタリングで利用可能なツール

ため、砂浜の保全・回復に多くの費用と時間を要するとして期待することができる。

1)

2)

ことが課題となっている。津波防災地域づくりと砂浜

国内ではCoastSatの適用事例は少ないため、その詳保全のあり方に関する懇談会（令和元年） では、後細についてはほとんど知られていない。そこで本技術

追い的な侵食対策から予測を重視した順応的な砂浜管資料では、Vosら 7)、8)およびCoastSatのマニュアル等理への展開を図ることが重要であると指摘しており、を参考にCoastSatの汀線抽出アルゴリズム実行環境のそのために必要な施策のひとつとして、砂浜幅のモニ構築および解析手順を説明するとともに、北海道沙流タリングを挙げている。

川河口で得られた測量データとの比較により、精度評

3)

砂浜幅を把握する際に一般的に用いられている方法

価行った結果を報告する。

はトータルステーションやレベル、GNSSなどの測量機器を用いた汀線測量である。汀線測量では、岸沖方

２．CoastSatの汀線抽出アルゴリズム

向に設定した測線上の高さと平面座標を複数点求める
ことにより砂浜の縦断地形が得られ、さらに、測線を２．１

衛星画像の前処理沿岸方向に複数設定することにより、汀線位置（例えCoastSatでは表－1に示すLandsat5、7、8、9および

ば、T.P. 0 mなど）を面的に把握することができる。

Sentinel-2の衛星画像から汀線を抽出する。各衛星画

汀線測量では汀線位置を高い精度で得ることができる

像はクラウドマスキングにより画素ごとに雲の有無がが、測量作業を人力で行う必要があるため、広範囲を判定され、雲の占める割合が閾値（４．２で説明する対象としたモニタリングでは多大な労力と費用を要す

cloud threshで指定）を超えた画像では汀線は抽出さることが課題となる。また、航空写真を用いた汀線のれない。

読み取りも砂浜幅の把握に有効といえる。航空写真か通常、衛星には、高解像度のモノクロ画像を撮影すらの汀線読み取りは広範囲の汀線位置を連続的に把握るパンクロマティックセンサーと低解像度のR（赤）、できるという点で測量よりも有利である。しかしながG（緑）、B（青）、NIR（近赤外線）、SWIR1（短波ら、汀線モニタリングを航空写真で行う場合、航空機赤外線）画像を撮影するマルチスペクトルセンサーによる写真撮影を定期的・継続的に行う必要があり、が備えられており、そのうち汀線抽出で使用する画多くの費用を要することが課題となる。

像は低解像度マルチスペクトルセンサー画像である。

近年、測量や航空写真からの読み取りよりも効率的
寒地土木研究所月報№851

2023年12月

マルチスペクトルセンサー画像の解像度はLandsat5、
49

表－1
衛星名

CoastSatで使用する衛星画像
運用期間

再訪周期

画像解像度

（年）

（日）

（m/pixel)

Landsat51984 - 2013

1630 (15)Landsat71999 - 現在

1630 (15)Landsat82013 - 現在

1630 (15)Landsat92021 - 現在

1630 (15)Sentinel-2

2015 - 現在

520 (10)※画像解像度のカッコ内はダウンサンプリング処理およびパンシャープン処理後の値である。
図－1

7、8、9で30 m、Sentinel-2で20 mである。これらの

Vosら7)のFig.4.(a)、(b)、(c)を掲載。オーストラ
リアNarrabeenの衛星画像（Sentinel-2）から得られた(a)RGB画像、(b)sand（オレンジ）、whitewater（水色）、water（青）、other land features（色付けなし）の階層に分類した結果、(c)MNDWI画像

低解像度画像に対してbi-linear補間によるダウンサンプリング処理や高解像度パンクロマティックセンサー画像を用いたパンシャープン処理を施すことにより、Landsat5、7、8、9とSentinel-2のマルチスペクトル画
像の解像度はそれぞれ15 mと10 mとなる。
各衛星の再訪周期はLandsat5、7、8、9で16日、Sentinel-2で5日であり、従来の汀線測量や航空写真の読み取りと比較して、かなり高い頻度で汀線を抽出す図－2ることが可能となる。
２．２

汀線位置の決定Vosら7)のFig.4.(e)を掲載。図－1(c)のMNDWI画像の確率密度のヒストグラム（オレンジ：sand、青：water、緑：whitewater、紫：other landfeatures、点線：大津の方法で得られたsandとwaterの境界

CoastSatでは、ニューラルネットワークを使用して、図－1に示すように、衛星画像の各画素をsand（砂浜）、water（海域）、whitewater（砕波による白波）、と重なるように分布する。汀線位置、すなわちsandother land features（その他陸上にあるもの）の4つのとwaterのMNDWIの境界は、汀線抽出精度の低下を階層に分類する。この階層分けで用いるパラメータは招くwhitewaterとother land featuresのMNDWI を除

オーストラリアの5つの海岸の衛星画像を学習データ

いたsandとwaterのみのMNDWIヒストグラムに対しとして事前に最適化されている 。て、画像の二値化のしきい値を求める際に広く適用さ

9)

4つの階層へ分類した後、各階層の画素のModified
Normalized Difference Water Index（修正正規化水指数：MNDWI）を計算する。MNDWIはSWIR1画像とG画像の強度を用いて以下の式で求められる。

れている大津の方法10)から得ることができる。
（図－1、図－2の点線）。
sandとwaterの境界値が得られた後、MNDWI画像に汀線として境界値の等値線を引く（図－1(c)参照）。
この作業ではマーチングスクウェア法11)が適用され、𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀 �𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆1 � �
𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆1 � �(1)

最終的にサブピクセルオーダー（Landsat5、7、8、9では15 m以下、Sentinel-2では10 m以下）の精度で汀線が引かれる。

MNDWIの確率密度をヒストグラムで描くと、図－2に示すように、sandとwaterのMNDWIはそれぞれ正２．３

汀線位置の潮位補正と負の領域にピークを有する二峰型の分布となる。ま

衛星画像から得られた汀線位置は図－3に示すようた、whitewaterのMNDWIは正から負の広い範囲にに、潮位のバイアスが含まれている。CoastSatでは基分布し、other land featuresのMNDWIはSandの範囲

準水位と撮影時の潮位との差および砂浜の勾配から、50

寒地土木研究所月報№851

2023年12月

て、解析パッケージをインストールする際に必要なソフトウェア（例えば、OpenSSL1.1 Full versionなど）が不足しているとエラーが表示される場合がある。その場合、Anaconda Promptでconda clean -allと入力して、途中までインストールされたパッケージのアンインストールを行い、不足しているソフトウェアをイン図－3

ストールした後、再度マニュアルに従って解析パッケ

汀線位置の潮位補正の概略図ージをインストールしなければならない。不足してい∆𝐿𝐿るソフトウェアが複数ある場合、この作業を繰り返し行う必要がある。また、CoastSatでは、Google EarthEngineを使用することから、事前にGoogleアカウン

以下の式により線形的に汀線位置を補正する。
∆𝐿𝐿𝑧𝑧���𝑧𝑧���

∆𝐿𝐿 �𝑧𝑧��� � 𝑧𝑧����𝑠𝑠 𝑧𝑧����

∆𝐿𝐿

トを用意しなければならない。
(2)

ここで、 ∆𝐿𝐿は汀線位置の補正量(m)、 𝑧𝑧��� は基準水位(m)、 𝑧𝑧����は撮影時水位(m)、 𝑠𝑠は砂浜の勾配を表す。

𝑧𝑧���
３．CoastSatの実行環境構築
𝑠𝑠𝑧𝑧����

CoasSatはPythonで記述されたオープンソースプ
𝑠𝑠𝑧𝑧����
ログラムであり、開発者のKilianVosによって最新CoastSatの実行環境を構築した後、Anaconda Prompt
でconda activate coastsatと入力して、Anacondaの環境（仮想環境）をBaseからcoastsatに変更する。その後、
earthengine authenticateと入力し、用意したGoogleアカウントに許可を与えればCoastSatを使用する準備が完了する。環境の変更とGoogleアカウントへの許可はCoastSatを使用する際に毎回必要となる。
CoastSatのメインプログラムの表示、修正、実行等
はPythonの統合開発環境（プログラムの記述、実行等を行うソフトウェア）であるJupyter NotebookやのプログラムがGitHub（URL: https://github.com/Spyderを使用する。Jupyter Notebook（Spyder）をkvos/CoastSat）で公開されている。CoastSatをイン
𝑠𝑠
ストールする際には、まず、Pythonのオープンソー使用する場合、Anaconda PromptでJupyter Notebook（Spyder）と入力してJupyter Notebook（Spyder）をスディストリビューションであるAnacondaをインス起動した後、メインプログラムを読み込み、適宜、修トールする必要がある。Anacondaのインストーラ正、実行等を行う。

ーはAnacondaのホームページ（URL: https://www.
anaconda.com）からダウンロードすることができる。

４．CoastSatの解析手順

Windows版の場合、インストーラーをダウンロードして実行すると自動的にAnacondaがインストールされCoastSatのメインプログラムはCoastSat-masterにると同時にPythonの実行環境が構築される。Python

保存されているサンプルプログラム（example.py）をの実行環境を構築したあと、GitHubからCoasSatのコ必要な個所のみ修正することによって作成することがード一式をZipファイル形式でダウンロードして、任できる。メインプログラムは解析画像の設定、汀線抽

意のフォルダで解凍するとCoastSat-masterというフ出パラメータの設定、測線の設定、潮位補正の設定のォルダが作成される。CoastSat-masterの中にはコー4つの項目から構成されている。以下に各項目の設定

ドやサンプルファイル、マニュアル（README.md

方法を説明する。

ファイル）などが保存されている。
CoastSat-masterが作成された後、Anaconda Prompt

４．１

解析画像の設定（または、Anaconda Powershell Prompt）を立ち上表－1に示すLandsat5、7、8、9およびSentinel-2のげて、cd（Change Directory）コマンドでCoastSat-中から使用する衛星画像を選択して（例えば、L5、masterに移動する。あとは、マニュアルに記載されL8、S2など複数選択可能）、期間と範囲（緯度、経度）ている通りにコマンドを入力して、必要な解析パッケを指定すると、指定した期間に撮影された全ての衛星ージをすべてインストールすればCoastSatの実行環境画像がGoogle Earth Engineから自動的にダウンローが構築される。なお、使用するパソコンの環境によっドされる。解析期間が長い場合、すべての衛星画像が

寒地土木研究所月報№851

2023年12月

51

ダウンロードされるまで半日以上を要する場合があ

４．３

測線の設定る。なお、マニュアルでは、解析範囲は100 km より測線は始点と終点の座標をメインプログラム上に直も小さくなるように設定することが推奨されており、接入力する、または、始点と終点の座標をgeojsonデ

解析範囲が100 km を超える場合は領域を分けて解析ータ形式で記載したファイルを読み込むことにより設することが望ましい。

定する。座標値は前項で説明したoutput_epsgの投影22座標系で入力する。なお、ArcGIS上で起動するソフ
４．２

トウェアDigital Shoreline Analysis System (DSAS)6)

汀線抽出パラメータの設定

汀線抽出パラメータは複数用意されており、そのうちユーザーによる設定が必要なパラメータはcloud_を使用することで、測線のgeojsonファイルを容易に作成することができる。

thresh、output_epsg、min_beach_area、min_length_sl、sand_color、max_dist_refである。cloud_threshは４．４

潮位補正の設定

ダウンロードした衛星画像に占める雲領域の割合のCoastSatでは、潮位データと砂浜の勾配から各測線許容値であり、この値を超える割合で雲が存在するの汀線位置を潮位補正した状態で出力することができ画像は解析されない。初期値は0.1と設定されている。

る（２．３参照）。潮位データは解析期間よりも長くとoutput_epsgは解析結果の出力で使用する投影座標系る必要があり、実測データまたはTPXO912)などの潮のEPSGコードである。例えば、解析結果をWGS84/位モデルからcsvファイルで作成する。砂浜の勾配のUTM zone 54Nで出力する場合、output_epsg は実測値が無い場合、Bujanら13)を参考に砂の粒径から32654となる。min_beach_area、min_length_slはそれ勾配を推定することができる。なお、CoastSatで入力ぞれ砂浜として認識する最小領域（単位：m ）と最できる勾配は1つであり、通常、平均的な勾配を入力小周長（単位：m）であり、良好な解析結果が得られする。勾配の時間変化を考慮する場合、CoastSatで潮るように、解析する領域の規模に応じて両値を調整位補正をせずに汀線位置を出力して、式(2)により個する。sand_colorは砂浜の色であり、default、latest、別に汀線位置を補正する必要がある。

2

dark、brightの中から選択して入力する。通常はdefaultを選択すればよいが、この設定で汀線が抽出

５．沙流川河口の汀線測量データを用いた精度評価されない場合、砂浜の色に応じて、dark、brightを選択する。latestはすべての砂浜の学習データから構築

５．１

汀線測量データ

されたパラメータであり、多くの砂浜に対して適応す汀線測量は図－4に示す沙流川河口の右岸と左岸のることができるが、汀線の抽出精度はあまり高くはな

砂浜において、2018年9月24日～2021年2月25日の間にい。

月1回程度の頻度で低潮位時に実施した（右岸：2018

CoastSatでは、汀線の誤抽出を低減するため、年9月24日～2021年2月22日の間に26回、左岸：2019年

Reference shoreline（参考汀線）と汀線抽出範囲を設

1月25日～2021年2月25日の間に24回）。測量ではGNSS定することができる。Reference shorelineを設定する

アンテナを持った状態で砂浜を3～4時間程度歩き、そ場合、雲が少なく汀線がはっきりと認識できる画像をの間、砂浜の位置、高さのデータを1秒間隔で記録し選択して、この画像上に目視で汀線位置を確認しながらその位置を手動（マウスクリック）で複数点プロットして設定する。また、汀線抽出範囲はReferenceshorelineからの距離をmax_dist_ref（単位：m）で設定する。Reference shorelineの設定はオプション機能であり、設定しなくても汀線を抽出することができるが、汀線の誤抽出を低減できるため利用することが望ましい。
ここで説明した6つのパラメータ以外にもユーザーが設定できるパラメータが用意されているが、それらに関しては基本的に初期値から変更する必要はない。
図－452

沙流川河口の位置図
寒地土木研究所月報№851

2023年12月

図－5

GNSS測量データから作成された50 cmメッシュのDEMの一例と測線（右岸測量日：2020年9月17日、左岸測量日：2020年9月16日）

図－7(a)右岸と(b)左岸の各測線における勾配の時間変化

表－2
右岸
測線名

図－6 形プロファイルからの勾配の求め方（図－5、測地線R5の地形プロファイル）た。測量終了後、kriging法を用いて、得られたデー

各測線における勾配の時間平均値
勾配

左岸
測線名

勾配

R10.094

L10.096

R20.094

L20.086

R30.092

L30.082

R40.099

L40.083

R50.095

L50.067

R60.094

L60.064

R70.092

L70.050

R80.089--R90.093--

５．４

砂浜の勾配データ

タから50 cmメッシュのDEM（T.P.基準）を作成し、５．１で作成したDEMから各測線における地形プロ

DEMから作成したT.P. = 0 mの等値線を汀線としたファイルを作成し、T.P. -0.5 m～T.P. +1.0 mの範囲で（図－5参照）。

勾配を求めた（図－6参照）。図－7と表－2はそれぞれ各測線における勾配の時系列と平均値を示している。

５．２

測線の設定

図－7では、右岸、左岸ともに、各測線における勾配測線は図－5に示すように、右岸ではR1～R9の9本、は時間的に変化している様子が確認できる。Vosら7)

左岸ではL1～L7の7本設定した。測線間隔はすべて50は勾配の時間変化を考慮して潮位補正を行った場合、mとした。

平均値を用いた場合と比べて、統計誤差に大きな改善が見られなかったことを示している。また、筆者も事

５．３

前に複数の測線において、勾配の時間変化を考慮して潮位データ

沙流川河口では潮位観測が行われていないため、潮潮位補正を行った結果と平均値を用いた結果で精度を位モデルTPXO9 を使用して、沙流川河口沖の水深比較し、両者に大きな違いがなかったことを確認して10 m地点における2018年8月1日～2021年3月1日の潮いる。したがって、本技術資料では、表－2に示す勾位を1時間間隔で求めた。

配の平均値を用いて汀線位置の潮位補正を行った結果

12)を報告する。
寒地土木研究所月報№851

2023年12月

53

図－8

Sentinel-2（撮影日：2020年12月8日）から得られた(a)RGB画像、(b)sand（オレンジ）、whitewater（水色）、water（青）、other land features（色付けなし）の階層に分類した結果、(c)MNDWI画像（砂浜幅が狭いため、(b)sand（オレンジ）が汀線に隠れている。）

図－9

 線R7における精度評価結果一覧（(a)汀線位置の時系列（青：測量データ、オレンジ：衛星画像）、(b)(a)のプロ
測ット（測量日の前後5日以内に撮影された衛星画像の結果、L7: Landsat7、L8:Landsat8、S2: Sentinel-2）、(c)汀
線位置の誤差、(d)誤差の確率密度

５．５

CoastSatの精度評価結果

図－9は測線R7における精度評価結果一覧であり、図－8は図－1と同様にSentinel-2から得られた(a)(a)汀線位置の時系列（青：測量データ、オレンジ：

RGB画像、(b)各画素をsand、whitewater、water、衛星画像）、(b)(a)のプロット、(c)汀線位置の誤差、other land featuresの階層に分類した画像、(c)(d)誤差の確率密度分布を示している。(b)のプロット

MNDWI画像を示している。なお、紙面の都合上、図では、測量日の前後5日以内に撮影された衛星画像のの掲載を割愛しているが、Sentinel-2よりも解像度が結果のみを示している。測線R7の精度は他の測線と低いLandsat8においても、図－8と同様に汀線が抽出比べてかなり高く、(a)に示す衛星画像から抽出されされていることを確認している。

た汀線位置は測量データから得られた汀線位置の変化

54

寒地土木研究所月報№851

2023年12月

図－10

測線R2における図－9と同様の図に対応しており、(b)に示すプロットでは相関係数が0.85を超えている。
図－10は測線R2における図－9と同様の図であり、他の測線と比べて精度が最も低かった結果を示している。測線R2では、測量データから得られた汀線は比較的変動が小さいのに対して、衛星画像から抽出された汀線は大きく変動しており、両者の相関係数は0.458とかなり低い値となる。
図－11は(a)右岸と(b)左岸の各測線における測量データと衛星画像から得られた汀線位置の相関係数（図－9と図－10の(b)に対応）を示している。右岸では、測線R1～R3の相関係数が比較的低い値を示す一方、測線が増えて、河川と海域の合流点に近づくにしたがって相関係数が増加する傾向が確認できる。左岸

図－11(a)右岸と(b)左岸の各測線における測量データと衛星画像から得られた汀線位置の相関係数

では、測線L3でやや低い値を示すが、それ以外の測線では相関係数に大きな違いは見られない。測線R1～R3の相関係数が低い原因は定かではないが、図－4、６．まとめ

図－5で確認できるように、河口右岸では、西側から河川と海域の合流点に近づくにしたがって砂浜の幅が本技術資料では、Vosら 7)、8)およびCoastSatのマニ増加している。測線R1～R3の付近ではとくに砂浜のュアル等を参考にCoastSatの汀線抽出アルゴリズム、幅が狭いため、潮位が高い状態で衛星画像が撮影され実行環境構築方法、解析手順を説明した。また、沙流た場合、砂浜の幅が画像の解像度（図－1参照）より川河口の測量データとの比較により、CoastSatの精度も狭くなり、汀線を正しく抽出できなかった可能性が評価を行った結果、測量データと衛星画像から抽出さある。

れた汀線位置の相関係数は右岸で0.46～0.87、左岸で0.61～0.81となった。砂浜の幅が狭い測線では比較的

寒地土木研究所月報№851

2023年12月

55

精度が低く、砂浜の幅が広がるにしたがって精度が高Baart, F., Donchyts, G., Aarninkhof, S.: The state

まる傾向が見られた。砂浜の幅が狭い場合、潮位が高of the world’s beaches, Sci Rep, pp.1-11, 2018.

い状態で撮影された衛星画像では砂浜が解像しきれず

7)に汀線の抽出精度が低下した可能性がある。

Vos. K., Harley. M.D., Splinter, K.D., Simmons,J.A., Turner, I.L.: Sub-annual to multi-decadal

現在、北海道沿岸の多くの砂浜海岸に対して

shoreline variability from publicly available

CoastSatを適用し、過去約40年間の汀線変化と波浪とsatellite imagery, Coast Eng, 150, pp. 160-174,の関係の評価を行っている。結果がまとまり次第、月2019.

報にて報告する予定である。

8)

Vos. K., Splinter, K.D., Harley. M.D., Simmons,J.A., Turner, I.L.: CoastSat: a Google Earthengine-enable Python toolkit to extract shorelines

参考文献

from publicly available satellite imagery, Environ
1)

一般財団法人沿岸技術研究センター：海岸保全施設の技術上の基準・同解説、公益社団法人日本港

Model Softw, 2019, 122, 104528, 2019.
9)

湾協会、2018.
2)

3)

4)

classification, Zenodo Digital Repository, 2019.

田中茂信、小荒井衛、深沢満：地形図の比較に

10) Otsu, N.: A threshold selection method from

よる全国の海岸線変化、海岸工学論文集、40巻、

gray-level histograms, IEEE TranSyst Man

pp.416-420、1993.

Cybern, 20, pp.62-66, 1979.

津波防災地域づくりと砂浜保全のあり方に関す

11) Cipolletti, M.P., Delrieux, C.A., Perillo, G.M.E.,

る懇談会：砂浜保全に関する中間とりまとめ、

Cintia Piccolo, M.: Superresolution border2019.

segmentation and measurement in remote

Kulei, T., Guneroglu, A., Karsli, F., Dihkan, M.:sensing images, Comput Geosci, 40, pp.87-96,

Automatic detection of shoreline change on

2012.

coastal Ramsar wetlans of Turkey, Ocean Eng.
5)

6)

Vos, K. Training data for: CoastSat image

12) Egbert, Gary D., and Svetlana Y. Erofeeva.:

38, pp.1141-1149, 2011.

Efficient inverse modeling of barotropic ocean

Hagenaars, G., de Vries, S., Luijendijk, A.P., detides, Journal of Atmospheric and Oceanic

Boer, W.P., Peniers, A.J.H.M.: On the accuracy

Technology, 19.2, pp.183-204, 2002.

of automated shoreline detection derived from13) Bujan, N., Cox, R.,Masselink, G.: From fine sandsatellite imagery: A case study of the sandto boulders: Examing the relationship between

motor mega-scale nourishment, Coast Eng. 133,

beach-face slope and sediment size, Marine

pp.113-125, 2018.

Geology, 417, 106012, 2019.

Luijendijk, A., Hagenaars, G., Ranasinghe, R.,

大塚 淳一
OTSUKA Junichi
寒地土木研究所
寒地水圏研究グループ
寒冷沿岸域チーム
主任研究員
博士（工学）

56

寒地土木研究所月報№851

2023年12月