勾配図を作成する

概要

f:id:Chiakikun:20211230000905p:plain

DEM

f:id:Chiakikun:20211230164842p:plain

勾配図

f:id:Chiakikun:20211230165105p:plain

勾配のヒストグラム

コード

using System;
using System.Runtime.InteropServices;

using DotSpatial.Data;
using DotSpatial.Data.Rasters.GdalExtension;

namespace Kobaizu
{
    class Program
    {
        [DllImport("kernel32.dll", CharSet = CharSet.Auto)]
        static extern bool SetDllDirectory(string lpPathName);

        static void Main(string[] args)
        {
            string dllpath = @"DotSpatialソース内のgdalが入っているフォルダ";
            string loadfilepath = @"DEM.tif";
            string savefilepath = @"保存先.tif";
            
            SetDllDirectory(dllpath);

            // GeoTiff読み込み
            GdalRasterProvider d = new GdalRasterProvider();
            IRaster src = d.Open(loadfilepath);

            int ncol = src.NumColumns;
            int nrow = src.NumRows;
            int band_num = src.NumBands;
            string prj = src.ProjectionString;
            double nodata = src.NoDataValue;

            double[] pGT = src.Bounds.AffineCoefficients;
            double xllcenter = pGT[0];
            double cellsize_x = pGT[1];
            double rotate1 = pGT[2];
            double yllcenter = pGT[3];
            double rotate2 = pGT[4];
            double cellsize_y = pGT[5];

            // GeoTiff書き出し先作成
            IRaster dst = Raster.CreateRaster(savefilepath, null, ncol, nrow, 1, typeof(float), new[] { string.Empty });
            dst.NoDataValue = nodata;
            dst.ProjectionString = prj;
            dst.Bounds = new RasterBounds(nrow, ncol, new double[] { xllcenter - cellsize_x / 2, cellsize_x, 0, yllcenter - cellsize_y / 2, 0, cellsize_y });

            Kobaizu(src.Value, nrow, ncol, dst.Value);

            dst.Save();

            Console.WriteLine("終了しました。");
            Console.ReadKey();
            return;
        }


        static void Kobaizu(IValueGrid src, int nrow, int ncol, IValueGrid dst)
        {
            double dx = 8.763;
            double dy = 12.348;

            for (int x = 0; x < ncol; x++)
            {
                for (int y = 0; y < nrow; y++)
                {
                    dst[y, x] = -9999;

                    if ((x < 1) || (y < 1) || (x >= ncol - 1) || (y >= nrow - 1)) continue;
                    double H11 = src[y - 1, x - 1]; if (H11 == -9999) continue;
                    double H12 = src[y - 1, x];     if (H12 == -9999) continue;
                    double H13 = src[y - 1, x + 1]; if (H13 == -9999) continue;
                    double H21 = src[y, x - 1];     if (H21 == -9999) continue;
                    double H23 = src[y, x + 1];     if (H23 == -9999) continue;
                    double H31 = src[y + 1, x - 1]; if (H31 == -9999) continue;
                    double H32 = src[y + 1, x];     if (H32 == -9999) continue;
                    double H33 = src[y + 1, x + 1]; if (H33 == -9999) continue;

                    double Ix = ((H11 + 2 * H21 + H31) - (H13 + 2 * H23 + H33)) / (8 * dx);
                    double Iy = ((H11 + 2 * H12 + H13) - (H31 + 2 * H32 + H33)) / (8 * dy);
                    double I = Math.Pow(Math.Atan(Math.Pow(Ix, 2) + Math.Pow(Iy, 2)), 0.5) * 180 / Math.PI;

                    dst[y, x] = I;
                }
            }
            return;
        }
    }
}

58、59行目の値は、参考資料のP3に載っていた縦方向、横方向の1グリッド分の距離です。

準備

f:id:Chiakikun:20211230001444p:plain

前回同様「DotSpatial.Data.dll」と「DotSpatial.Data.Rasters.GdalExtension.dll」を参照に追加してからソースをコンパイルしてください。

出力

下図は勾配が50~60°の箇所を赤く表示させた勾配図に地理院地図を重ねてみた図です。google ストリートビューで黒丸の辺りから北を見てみると、赤い個所の急勾配具合が良くわかりました。

勾配図に地理院地図を重ねた 出典:国土地理院 (https://cyberjapandata.gsi.go.jp

 

最後までご覧いただき、ありがとうございました。

 

 





傾斜量図を作成する

概要

f:id:Chiakikun:20211230000905p:plain

DEM

f:id:Chiakikun:20211229232031p:plain

傾斜量図

コード

using System;
using System.Runtime.InteropServices;

using DotSpatial.Data;
using DotSpatial.Data.Rasters.GdalExtension;

namespace Keisharyozu
{
    class Program
    {
        [DllImport("kernel32.dll", CharSet = CharSet.Auto)]
        static extern bool SetDllDirectory(string lpPathName);

        static void Main(string[] args)
        {
            string dllpath = @"DotSpatialソース内のgdalが入っているフォルダ";
            string loadfilepath = @"DEM.tif";
            string savefilepath = @"保存先.tif";

            SetDllDirectory(dllpath);

            // GeoTiff読み込み
            GdalRasterProvider d = new GdalRasterProvider();
            IRaster src = d.Open(loadfilepath);

            int ncol = src.NumColumns;
            int nrow = src.NumRows;
            int band_num = src.NumBands;
            string prj = src.ProjectionString;
            double nodata = src.NoDataValue;

            double[] pGT = src.Bounds.AffineCoefficients;
            double xllcenter = pGT[0];
            double cellsize_x = pGT[1];
            double rotate1 = pGT[2];
            double yllcenter = pGT[3];
            double rotate2 = pGT[4];
            double cellsize_y = pGT[5];

            // GeoTiff書き出し先作成
            IRaster dst = Raster.CreateRaster(savefilepath, null, ncol, nrow, 1, typeof(float), new[] { string.Empty });
            dst.NoDataValue = nodata;
            dst.ProjectionString = prj;
            dst.Bounds = new RasterBounds(nrow, ncol, new double[] { xllcenter - cellsize_x / 2, cellsize_x, 0, yllcenter - cellsize_y / 2, 0, cellsize_y });

            Keisharyozu(src.Value, nrow, ncol, dst.Value);

            dst.Save();

            Console.WriteLine("終了しました。");
            Console.ReadKey();
            return;
        }


        static void Keisharyozu(IValueGrid src, int nrow, int ncol, IValueGrid dst)
        {
            double dx = 8.763;
            double dy = 12.348;

            for (int x = 0; x < ncol; x++)
            {
                for (int y = 0; y < nrow; y++)
                {
                    dst[y, x] = -9999;

                    if ((x < 1) || (y < 1) || (x >= ncol - 1) || (y >= nrow - 1)) continue;
                    double H11 = src[y - 1, x - 1]; if (H11 == -9999) continue;
                    double H12 = src[y - 1, x];     if (H12 == -9999) continue;
                    double H13 = src[y - 1, x + 1]; if (H13 == -9999) continue;
                    double H21 = src[y, x - 1];     if (H21 == -9999) continue;
                    double H23 = src[y, x + 1];     if (H23 == -9999) continue;
                    double H31 = src[y + 1, x - 1]; if (H31 == -9999) continue;
                    double H32 = src[y + 1, x];     if (H32 == -9999) continue;
                    double H33 = src[y + 1, x + 1]; if (H33 == -9999) continue;

                    double Sx = (H11 + H21 + H31 - (H13 + H23 + H33)) / (6 * dx);
                    double Sy = (H11 + H12 + H13 - (H31 + H32 + H33)) / (6 * dy);
                    double S = Math.Pow((Math.Pow(Sx, 2) + Math.Pow(Sy, 2)), 0.5);

                    dst[y, x] = S * 1000;
                }
            }
            return;
        }
    }
}

58、59行目の値は、参考資料のP3に載っていた縦方向、横方向の1グリッド分の距離です。

準備

f:id:Chiakikun:20211230001444p:plain

VisualStudioの参照の追加で「DotSpatial.Data.dll」と「DotSpatial.Data.Rasters.GdalExtension.dll」を追加してからソースをコンパイルします。

出力結果

下図は傾斜量図と、地理院地図を重ねた図になります。赤が濃く出ている箇所ほど傾斜が強くなっています。ストリートビューで「幌延町 桂橋」と検索して移動した地点で北に向くと、赤い個所の地形がどんなものかよくわかると思います。。

傾斜量図

地理院地図を重ねた 出典:国土地理院 (https://cyberjapandata.gsi.go.jp

 

最後までご覧いただき、ありがとうございました。

 

 





QGISのプラグイン作成 PostGISのテーブルをマップに表示する

はじめに

PyQGISでPostGISのテーブルをマップに追加する方法について書いています。
このページでは国土数値情報から、避難所と行政区域をPostGISに取り込んで使っています。
行政区域は県毎にレコードを結合しました。


テーブル丸ごと開く

    
        uri = qgis.core.QgsDataSourceUri()
        uri.setConnection(サーバーアドレス, ポート, データベース名, ユーザーID, パスワード)

        uri.setDataSource("", "避難所", "geom")
        vlayer = QgsVectorLayer(uri.uri(),"避難所","postgres")
        if not vlayer.isValid():
            QMessageBox.about(None, "エラー", "開けませんでした。")
        else:
            QgsProject.instance().addMapLayer(vlayer)

f:id:Chiakikun:20210521225819p:plain
shp等のファイルと変わらず、マップに追加できます。


クエリを開く

行政区域_都道府県(国土数値情報の行政区域のオブジェクトを県毎に結合して加工したレイヤです)の千葉県に交差する避難所をレイヤとして追加します。

    
        uri = qgis.core.QgsDataSourceUri()
        uri.setConnection(サーバーアドレス, ポート, データベース名, ユーザーID, パスワード)

        sql = "select 避難所.* from 避難所, 行政区域_都道府県 where 行政区域_都道府県."N03_001" = '千葉県' and st_intersects(避難所.geom, 行政区域_都道府県.geom)"
        uri.setDataSource("",f'({sql})',"geom", "", "id")
        vlayer = QgsVectorLayer(uri.uri(),"避難所","postgres")
        if not vlayer.isValid():
            QMessageBox.about(None, "エラー", "開けませんでした。")
        else:
            QgsProject.instance().addMapLayer(vlayer)

f:id:Chiakikun:20210521230931p:plain

このように、特定のレコードだけ抽出して開くことができます。ただ...

 

f:id:Chiakikun:20210521231051p:plain

編集はできないようです。


Viewを開く

まずPostgreSQLで、千葉県に交差する避難所を抽出するViewを作成します。

CREATE VIEW 避難所ビュー AS
SELECT
    避難所.*
FROM
    避難所,
    行政区域_都道府県
WHERE
    行政区域_都道府県."N03_001" = '千葉県'
AND
    st_intersects(避難所.geom, 行政区域_都道府県.geom);

Viewを作成できたら、実行します。

    
        uri = qgis.core.QgsDataSourceUri()
        uri.setConnection(サーバーアドレス, ポート, データベース名, ユーザーID, パスワード)

        uri.setDataSource('public', '避難所ビュー', 'geom', '', 'id')
        vlayer = QgsVectorLayer(uri.uri(False),"避難所ビュー","postgres")
        if not vlayer.isValid():
            QMessageBox.about(None, "エラー", "開けませんでした。")

f:id:Chiakikun:20210521232023p:plain
この方法も千葉県に交差する避難所のレコードだけのレイヤが追加され...

 

f:id:Chiakikun:20210521232142p:plain

一見編集も可能に見えますが...

 

f:id:Chiakikun:20210522010101p:plain

編集を確定しようとすると「ジオメトリの値変更中のPostGISエラー」とエラー表示されてしまいます。編集可能なViewは、単一のテーブルから作らないといけないようなので、千葉県内の避難所をViewで抽出するには...

    
CREATE VIEW 避難所ビュー AS
SELECT
    避難所.*
FROM
    避難所
WHERE
    LEFT("p20_003", 3) = '千葉県';

避難所の「p20_003」フィールドには住所が入っているので、この値を使ってViewを作成します。このVIewであれば編集が可能です。


使い道

クエリやViewを使えば、手間をかけて県毎にわける必要が無くなり、プラグインを作りこめばボタン一発で必要箇所のレコードだけをレイヤとして追加できるようになります。
また、Viewなら複数人で一つのテーブルの編集作業をする際に、Aさんは千葉県、Bさんは埼玉県...という感じで分けられるので、同じレコードを編集してしまうことも無くなるのではと思います(まだ実際に複数人でビューの編集を試していないので不安ですが)。


最後までご覧いただき、ありがとうございました。