コンテンツにスキップ

FAQ/ 占有率のラスターデータの作成

あるメッシュに占める田んぼの占有面積を求めたいです。 方法を教えてください。

回答

DioVISTAにそうした機能はありません。

ご要望を、QGISで実現する方法を示します。

ここでは、田んぼの占有率を求める例を示します。同じ方法は、建物占有率を求める場合にも応用できます。

田んぼのポリゴンデータとして、農林水産省 筆ポリゴン 公開サイト からダウンロードしたデータを使います。このデータは、GeoJson形式のポリゴンです。

  1. QGISを立ち上げます。

    1. メニュー > [プロジェクト] > [新規] を選択します。
  2. 田んぼポリゴンを読み込みます。

    1. メニュー > [レイヤ] > [レイヤを追加] > [ベクタレイヤを追加] を選択します。
    2. ダイアログが出現します。
    3. [ベクタデータセット] に、田んぼのポリゴンファイルを指定します。ここでは、 2022_151017.json を指定しました。
    4. [追加]ボタンを押します。
    5. [閉じる]ボタンを押します。
    6. 田んぼポリゴンが読み込まれました。
  3. メッシュのポリゴンを作ります。ここでは、100mメッシュを作ります。

    1. メニュー > [ベクタ] > [調査ツール] > [グリッドを作成] を選択します。
    2. ダイアログが出現します。
    3. [グリッドタイプ] を [長方形(Polygon)] にします。
    4. [グリッドの範囲] を指定します。これは、3次メッシュと一致するようにします。経度方向は 1/80の倍数、緯度方向は 1/120の倍数にします。この例では、[139.125,139.375,37.75,38] にしました。
    5. [水平方向の間隔] を [0.001250] にします。この値の求め方: [メッシュサイズ]/80000
    6. [鉛直方向の間隔] を [0.000833]にします。この値の求め方: [メッシュサイズ]/120000
    7. [実行]を押します。
    8. [閉じる]を押します。
    9. グリッドポリゴンが作成されました。
  4. 田んぼポリゴンをメッシュで分割します。

    1. メニュー > [プロセッシング] > [ツールボックス] で、画面左側に プロセッシングツールボックスを表示します。
    2. 画面左側 プロセッシングツールボックス > [ベクタオーバーレイ] > [交差 (intersect)] を選択します。
    3. ダイアログが出現します。
    4. [入力レイヤ] を 田んぼポリゴンのレイヤにします。この例では、[2022_151017] にしました。
    5. [オーバーレイレイヤ] を [グリッド] にします。
    6. [実行]を押します。
    7. [閉じる]を押します。
    8. 交差レイヤーが作成されました。
  5. 分割した田んぼポリゴンをメッシュごとに統合します。

    1. 画面左側 プロセッシングツールボックス > [ベクタジオメトリ] > [融合 (dissolve)] を選択します。
    2. ダイアログが出現します。
    3. [入力レイヤ] を、[交差(intersect)] にします。
    4. [基準となる属性]を、[id]にします。
    5. [実行]を押します。
    6. [閉じる]を押します。
  6. メッシュごとの田んぼの面積を求めます。

    1. レイヤの [融合ポリゴンの出力] を右クリックし、[属性テーブルを開く] を選択します。
    2. フィールド計算機を開くボタン を押します。
    3. ダイアログが出現します。
    4. [出力する属性(フィールド)の名前] を [t_area] にします。
    5. [フィールド型] を [小数点付き数値(real)] にします。
    6. [フィールド長] を [12]に、[精度] を [6] にします。
    7. [式] に [$area] と入力します。
    8. [OK]を押します。
    9. 属性テーブルを閉じます。
  7. メッシュの属性に、田んぼの面積を追加します。

    1. 画面左側 プロセッシングツールボックス > [ベクタ一般] > [属性テーブルで結合 (table join)] を選択します。
    2. ダイアログが出現します。
    3. [入力レイヤ] を [グリッド] にします。
    4. [入力レイヤの結合対象フィールド] を [id] にします。
    5. [第2の入力レイヤ] を [融合ポリゴンの出力] にします。
    6. [第2の入力レイヤの結合対象フィールド] を [id] にします。
    7. [第2の入力レイヤからコピーする属性] を [t_area] にします。
    8. [実行]を押します。
    9. ログに赤文字で下記のようなメッセージが出ることがあります。

      入力レイヤの 51657 個の地物が一致しませんでした

      これは、いくつかのメッシュでは、その中に1つも田んぼポリゴンがないことを伝えています。これは問題ないので無視します。 1. [閉じる]を押します。

  8. メッシュの属性に、田んぼ占有率を追加します。

    1. レイヤの [出力レイヤ] を右クリックし、[属性テーブルを開く] を選択します。
    2. フィールド計算機を開くボタン を押します。
    3. ダイアログが出現します。
    4. [出力する属性(フィールド)の名前] を [t_rate] にします。
    5. [フィールド型] を [小数点付き数値(real)] にします。
    6. [フィールド長] を [12]に、[精度] を [6] にします。
    7. [式] に下記の式を記載します。

      case when "t_area" is null then 0 else "t_area"/$area end

    8. [OK]を押します。

    9. 属性テーブルを閉じます。
  9. 作成したデータを可視化します。

    1. レイヤの[出力レイヤ]を右クリックし、[プロパティ]を選択します。
    2. ダイアログが出現します。
    3. 左側のタブから [シンボロジ] を選択します。
    4. [連続値による定義(graduated)] を選択します。
    5. [値] を [t_rate] にします。
    6. [シンボル] をクリックし [不透明度] を [50.0%] にします。
    7. [分類]をクリックします。白から赤に変化する5種類の色が定義されました。
    8. [OK]を押します。
    9. 可視化された田んぼ占有率を示します。
  10. 作成したデータをラスタデータに変換します。

    1. メニュー > [ラスタ] > [変換] > [ベクタをラスタ化(rasterlize)] を選択します。
    2. ダイアログが出現します。
    3. [入力レイヤ] を [出力レイヤ] にします。
    4. [焼き込む値の属性(フィールド)] を [t_rate] にします。
    5. [出力ラスタサイズの単位] を [地理単位] にします。
    6. [水平方向の解像度] を ステップ3の[水平方向の間隔]と一致させます。この例では、[0.001250] にしました。
    7. [鉛直方向の解像度] を ステップ3の[鉛直方向の間隔]と一致させます。この例では、[0.000833] にしました。
    8. [出力領域] ステップ3の[グリッドの範囲]と一致させます。この例では、[139.125,139.375,37.75,38] にしました。
    9. [実行]を押します。
    10. [閉じる]を押します。
  11. 作成したラスタデータをESRI Ascii Raster形式で保存します。

    1. メニュー > [ラスタ] > [変換] > [形式変換(gdal_translate)] を選択します。
    2. ダイアログが出現します。
    3. [入力レイヤ] を [出力ファイル] にします。
    4. [出力ファイル] を 拡張子が *.asc のファイルとします。ここでは、[test_t_rate.asc] としました。
    5. [実行]を押します。
    6. ログに 下記のようなメッセージが出ることがあります。

      Warning 1: Producing a Golden Surfer style file with DX and DY instead of CELLSIZE since the input pixels are non-square. Use the FORCE_CELLSIZE=TRUE creation option to force use of DX for even though this will be distorted. Most ASCII Grid readers (ArcGIS included) do not support the DX and DY parameters.

      これは、緯度と経度の解像度が異なるため、本来の ESRI Ascii Raster 形式では表現できず、DX, DYというパラメータを使うことを伝えています。これは問題ないので無視します。

    7. [閉じる]を押します。

    8. 作成されたファイル [test_t_rate.asc] の先頭から7行を示します。
      1
      2
      3
      4
      5
      6
      7
      ncols        200
      nrows        300
      xllcorner    139.125000000000
      yllcorner    37.750000000000
      dx           0.001250000000
      dy           0.000833333333
      NODATA_value  0
      

関連項目


最終更新日: 2023-07-12