とらりもん - 層別ランダムサンプリング Diff
- Added parts are displayed like this.
- Deleted parts are displayed
like this.
{{toc}}
! 理論
[[データとデータ解析,2.層別サンプリング|http://ebsa.ism.ac.jp/ebooks/sites/default/files/ebook/1824/pdf/ch08-02.pdf]]
! 実践
日本域土地利用・土地被覆図を例に，試す．
# first, merge all .tif file
$ gdal_merge.py v16.09/*.tif [ここでtab] -o Japan.tif
!! 各カテゴリのサンプル数を指定する方法（採用）
GRASSの[[r.sample.category|https://grass.osgeo.org/grass72/manuals/addons/r.sample.category.html]]で良さそう．
* 内部では[[r.random|https://trac.osgeo.org/grass/browser/grass/branches/releasebranch_7_2/raster/r.random]]によって点を発生させている．
* r.randomでは[[G_lrand48|https://grass.osgeo.org/programming7/lrand48_8c.html]]が呼ばれている.
$ g.extension extension=r.sample.category operation=add # install r.sample.category
$ r.in.gdal input=Japan.tif output=Japan # Japan.tifが103GBなので，めちゃくちゃ時間かかる
$ g.region rast=Japan
$ sample_number=2000
$ date && r.sample.category input=Japan output=Japan_2000points npoints=2000 && date # めちゃくちゃ時間がかかる
$ v.out.ogr input=Japan_2000points format=KML output=Japan_2000points.kml # save as kml
参考までに，sample_number=2000の場合，かかった時間．
*r.sample.category：15時間30分
!! Google Earthの最新画像のみを用いる場合
* [[Laco-Wiki|https://laco-wiki.net/en/Welcome]]
!! 各カテゴリのサンプル数を指定せず，適当にばらまく方法（ボツ）
手順は以下の通り．
各カテゴリのサンプル数をm個とする．
# 日本域の緯度経度の範囲内に，乱数をn個（n >> m）発生させる（乱数にはインデックスを付けておく）.
# gdallocationinfoを用いて，カテゴリをインデックス順に各m個抽出する．
# 抽出した点をkmlに変換する．
# 変換したkmlをGoogle Earthにインポートして，目視チェック．
# 全ての点数をチェックし終えた後に，混合行列を作成する．
参考までに，m=1億の時の計測時間．
*乱数1億個発生：5分
*gdallocationinfo：4000点で3分なので，1億点で75000分（＝1250時間＝52日)
**ただし直列実行時．
! 混同行列の作成
$ mkdir -p shp kml_category # 必要なディレクトリを作成
GRASS GISで作成したKMLはサイズが大きいので，以下の変換を行う．
# バッファを作成する．
#カテゴリ0を取り除く
#カテゴリごとにファイルを分割する
QGIS上で行う．
各属性値は以下の通り．
*cat:id
*japan:カテゴリー（0-10）
!!! ファイルの読み込み
レイヤ -> ベクタ -> ベクタ読み込み
!!! 距離バッファを作成
今回検証する日本域土地被覆図の解像度は10 mなので，半径5 mに設定する．
#ベクタ -> 空間演算ツール -> 固定距離バッファ
距離：4.545454545454545e-05 # 1deg = 110kmとした
ファイル名：buffer
RUN
============================
CRSをメートルに変換してからバッファを計算したかったが，以下の手順ではなぜかメートルにならず．
#距離の単位をmに変更する
設定 -> カスタムCRS -> + -> 名称：WG84_meter -> パラメーターに以下を入力
+proj=longlat +datum=WGS84 +no_defs +units=m
OK.
*レイヤパネルで右クリック -> レイヤのCRSの設定 -> WGS84_meterを選択 -> OK
*レイヤパネルで右クリック -> レイヤのCRSをプロジェクトのCRSに設定
距離：5
RUN
WGS84_meter
OK
============================
!!! カテゴリごとにファイルを分割
計算結果は自動的にQGISに読み込まれるが，そのままでは文字コードの関係？でベクタレイヤの分割が行えない．
そこで既に読みこまれている「バッファ」を削除し，先ほど保存した「buffer」を再度読み込ませる．
ベクタ->データマネジメントツール->ベクタレイヤの分割
*ユニークIDフィールド：japan
*出力ディレクトリ：shp
Run
選択した出力ディレクトリに，buffer_japan_?.shp というシェープファイルが作成される．
!!! kmlに再変換
以下のスクリプトをshp2kml.shという名前で保存し，kmlに変換．
#!/bin/sh
shp_dir="shp"
kml_dir="kml_category"
for i in `seq 1 10`; do
ogr2ogr -f KML ${kml_dir}/category_${i}.kml ${shp_dir}/buffer_japan_${i}.shp
done
$./shp2kml.sh
!!Google Earthによる確認
作成したKMLをGoogle Earth上にインポートし，実際に精度検証を行う．
ファイル -> 開く -> Open.
* 過去画像への切り替え
表示 -> 過去のイメージ
* ID順に飛んでいく
! 理論
[[データとデータ解析,2.層別サンプリング|http://ebsa.ism.ac.jp/ebooks/sites/default/files/ebook/1824/pdf/ch08-02.pdf]]
! 実践
日本域土地利用・土地被覆図を例に，試す．
# first, merge all .tif file
$ gdal_merge.py v16.09/*.tif [ここでtab] -o Japan.tif
!! 各カテゴリのサンプル数を指定する方法（採用）
GRASSの[[r.sample.category|https://grass.osgeo.org/grass72/manuals/addons/r.sample.category.html]]で良さそう．
* 内部では[[r.random|https://trac.osgeo.org/grass/browser/grass/branches/releasebranch_7_2/raster/r.random]]によって点を発生させている．
* r.randomでは[[G_lrand48|https://grass.osgeo.org/programming7/lrand48_8c.html]]が呼ばれている.
$ g.extension extension=r.sample.category operation=add # install r.sample.category
$ r.in.gdal input=Japan.tif output=Japan # Japan.tifが103GBなので，めちゃくちゃ時間かかる
$ g.region rast=Japan
$ sample_number=2000
$ date && r.sample.category input=Japan output=Japan_2000points npoints=2000 && date # めちゃくちゃ時間がかかる
$ v.out.ogr input=Japan_2000points format=KML output=Japan_2000points.kml # save as kml
参考までに，sample_number=2000の場合，かかった時間．
*r.sample.category：15時間30分
!! Google Earthの最新画像のみを用いる場合
* [[Laco-Wiki|https://laco-wiki.net/en/Welcome]]
!! 各カテゴリのサンプル数を指定せず，適当にばらまく方法（ボツ）
手順は以下の通り．
各カテゴリのサンプル数をm個とする．
# 日本域の緯度経度の範囲内に，乱数をn個（n >> m）発生させる（乱数にはインデックスを付けておく）.
# gdallocationinfoを用いて，カテゴリをインデックス順に各m個抽出する．
# 抽出した点をkmlに変換する．
# 変換したkmlをGoogle Earthにインポートして，目視チェック．
# 全ての点数をチェックし終えた後に，混合行列を作成する．
参考までに，m=1億の時の計測時間．
*乱数1億個発生：5分
*gdallocationinfo：4000点で3分なので，1億点で75000分（＝1250時間＝52日)
**ただし直列実行時．
! 混同行列の作成
$ mkdir -p shp kml_category # 必要なディレクトリを作成
GRASS GISで作成したKMLはサイズが大きいので，以下の変換を行う．
# バッファを作成する．
#カテゴリ0を取り除く
#カテゴリごとにファイルを分割する
QGIS上で行う．
各属性値は以下の通り．
*cat:id
*japan:カテゴリー（0-10）
!!! ファイルの読み込み
レイヤ -> ベクタ -> ベクタ読み込み
!!! 距離バッファを作成
今回検証する日本域土地被覆図の解像度は10 mなので，半径5 mに設定する．
#ベクタ -> 空間演算ツール -> 固定距離バッファ
距離：4.545454545454545e-05 # 1deg = 110kmとした
ファイル名：buffer
RUN
============================
CRSをメートルに変換してからバッファを計算したかったが，以下の手順ではなぜかメートルにならず．
#距離の単位をmに変更する
設定 -> カスタムCRS -> + -> 名称：WG84_meter -> パラメーターに以下を入力
+proj=longlat +datum=WGS84 +no_defs +units=m
OK.
*レイヤパネルで右クリック -> レイヤのCRSの設定 -> WGS84_meterを選択 -> OK
*レイヤパネルで右クリック -> レイヤのCRSをプロジェクトのCRSに設定
距離：5
RUN
WGS84_meter
OK
============================
!!! カテゴリごとにファイルを分割
計算結果は自動的にQGISに読み込まれるが，そのままでは文字コードの関係？でベクタレイヤの分割が行えない．
そこで既に読みこまれている「バッファ」を削除し，先ほど保存した「buffer」を再度読み込ませる．
ベクタ->データマネジメントツール->ベクタレイヤの分割
*ユニークIDフィールド：japan
*出力ディレクトリ：shp
Run
選択した出力ディレクトリに，buffer_japan_?.shp というシェープファイルが作成される．
!!! kmlに再変換
以下のスクリプトをshp2kml.shという名前で保存し，kmlに変換．
#!/bin/sh
shp_dir="shp"
kml_dir="kml_category"
for i in `seq 1 10`; do
ogr2ogr -f KML ${kml_dir}/category_${i}.kml ${shp_dir}/buffer_japan_${i}.shp
done
$./shp2kml.sh
!!Google Earthによる確認
作成したKMLをGoogle Earth上にインポートし，実際に精度検証を行う．
ファイル -> 開く -> Open.
* 過去画像への切り替え
表示 -> 過去のイメージ
* ID順に飛んでいく