2010年12月9日木曜日

GMTでポリゴンデータを用いてグリッドデータをマスキング

GMT (The Generic Mapping Tools)でポリゴンデータで、グリッドデータをマスキングのやりかた。

以下の処理はDEMで計算した標高変化の分布図から、関心のある領域(例えば特定の氷河のポリゴンなど)だけを取り出す場合。

変数の説明
【グリッドデータ】
${trend_nc}:標高変化の分布データ
${trend_nc_cut}:標高変化の分布データを処理対象領域のみを切り抜いたデータ
${mask}:マスクグリッド(1か0の2値データ:ポリゴン内部が1、外部は0)
${mask_na}:マスクグリッドの0をNAN値に置き換えたもの
${trend_nc_masked}:標高変化の分布データのうち、ポリゴンの内部のみを切り抜いたデータ(←一番欲しいデータ)

【ベクター】
${GLACIER_OUTLINE}:くり抜きに用いるポリゴンデータ(GMT形式)

【その他】
${utm_range}:処理対象範囲(今回の例ではUTM座標)
${image}:GMTの出力画像ファイル
${trend_cpt}:makecptなどで生成したカラーパレットファイル


#標高変化のデータから処理対象流域のみを切り抜く
#(マスクグリッドと範囲を同一にするため)
grdcut ${trend_nc} -R${utm_range} -G${trend_nc_cut}

#ポリゴンからマスクグリッドの作成
grdmask -F -m ${GLACIER_OUTLINE} -G${mask} -R${utm_range} -I${CELLSIZE}

#マスクグリッドの0の値をNANの値に変更
grdmath ${mask}  0 NAN = ${mask_na}

#標高変化のデータ×マスクグリッドを計算することによって、
#ポリゴンの外部の値を全てNANに変える。
grdmath -N ${mask_na} ${trend_nc_cut} MUL = ${trend_nc_masked}

#出来上がったデータから図の出力
grdimage ${trend_nc_masked} -C${trend_cpt} -Q -R -J -K > ${image}

0 件のコメント:

コメントを投稿