DotSpatialを使って遷急線・遷緩線図と尾根・谷線図も作成してみました。

概要

  • 前回に引き続き、DotSpatialでDEMを読み込んで遷急線・遷緩線図等を作成します。
  • 参考資料に載っていた図の作成方法は少し大変そうだったので、大分簡略しました。

f:id:Chiakikun:20220101143723p:plain

DEM

f:id:Chiakikun:20220101143913p:plain

遷急線・遷緩線図

f:id:Chiakikun:20220101144146p:plain

尾根・谷線図

コード

using System;
using System.Runtime.InteropServices;
using System.Linq;

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

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

        static void Main(string[] args)
        {
            SetDllDirectory("DotSpatial.Data.Rasters.GdalExtensionのフォルダ");

            string loadfilepath = "DEMファイル.tif";

            // 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];

            // 遷急線・遷緩線図
            string savefilepath = "出力ファイル1.tif";
            IRaster kankyu = Raster.CreateRaster(savefilepath, null, ncol, nrow, 1, typeof(float), new[] { string.Empty });
            kankyu.NoDataValue = nodata;
            kankyu.ProjectionString = prj;
            kankyu.Bounds = new RasterBounds(nrow, ncol, new double[] { xllcenter - cellsize_x / 2, cellsize_x, 0, yllcenter - cellsize_x / 2, 0, cellsize_y });
            // 尾根・谷線図
            savefilepath = "出力ファイル2.tif";
            IRaster onetani = Raster.CreateRaster(savefilepath, null, ncol, nrow, 1, typeof(float), new[] { string.Empty });
            onetani.NoDataValue = nodata;
            onetani.ProjectionString = prj;
            onetani.Bounds = new RasterBounds(nrow, ncol, new double[] { xllcenter - cellsize_x / 2, cellsize_x, 0, yllcenter - cellsize_x / 2, 0, cellsize_y });

            Keishahenkansenzu(src.Value, nrow - 1, ncol - 1, kankyu.Value, onetani.Value);
            kankyu.Save();
            onetani.Save();

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


        static void Keishahenkansenzu(IValueGrid src, int nrow, int ncol, IValueGrid kankyu, IValueGrid onetani)
        {
            double dx = 8.763;
            double dy = 12.348;
            double dx2 = Math.Pow(dx, 2);
            double dy2 = Math.Pow(dy, 2);

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

            for (int x = 1; x < ncol; x++)
            {
                for (int y = 1; y < nrow; y++)
                {
                    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 H22 = src[y, x];         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 tate = (H12 - 2 * H22 + H32) / dx2;              // |
                    double yoko = (H21 - 2 * H22 + H23) / dy2;              // ー
                    double naname1 = (H11 - 2 * H22 + H33) / (dx2 + dy2);   // \
                    double naname2 = (H13 - 2 * H22 + H31) / (dx2 + dy2);   // /

                    double[] array_kyokuritsu = new double[] { Math.Abs(tate), Math.Abs(naname1), Math.Abs(yoko), Math.Abs(naname2) };

                    switch(array_kyokuritsu.ToList().IndexOf(array_kyokuritsu.Max()))
                    {
                        case 0:
                            if (((H12 <= H22 && H32 <= H22) || (H12 >= H22 && H32 >= H22)))
                                onetani[y, x] = tate;
                            else
                                kankyu[y, x] = tate;
                            break;
                        case 1:
                            if (((H11 <= H22 && H33 <= H22) || (H11 >= H22 && H33 >= H22)))
                                onetani[y, x] = naname1;
                            else
                                kankyu[y, x] = naname1;
                            break;
                        case 2:
                            if (((H21 <= H22 && H23 <= H22) || (H21 >= H22 && H23 >= H22)))
                                onetani[y, x] = yoko;
                            else
                                kankyu[y, x] = yoko;
                            break;
                        case 3:
                            if (((H13 <= H22 && H31 <= H22) || (H13 >= H22 && H31 >= H22)))
                                onetani[y, x] = naname2;
                            else
                                kankyu[y, x] = naname2;
                            break;
                    }

                    // この辺は好みに応じて調整する
                    if (-0.01 < onetani[y, x] && onetani[y, x] < 0.01) onetani[y, x] = -9999;
                    if (-0.01 < kankyu[y, x] && kankyu[y, x] < 0.01) kankyu[y, x] = -9999;                    
                }
            }
            return;
        }
    }
}

準備

f:id:Chiakikun:20211230001444p:plain

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

 

 

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

DotSpatialを使って曲率図も作成してみました。

概要

  • 前回に続いて、DEMから曲率図を作成します。
  • 今回もラスタの読み書きにはDotSpatialを使用しました。
  • 今月作成したページは、画像とコードを入れ替えてるだけの手抜き仕様。

f:id:Chiakikun:20211230000905p:plain

DEM

f:id:Chiakikun:20211231001109p:plain

曲率図

コード

using System;
using System.Runtime.InteropServices;

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

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

        static void Main(string[] args)
        {
            SetDllDirectory("DotSpatial.Data.Rasters.GdalExtensionのパス");

            string loadfilepath = "DEM.tif";
            string savefilepath = "出力ファイル.tif";

            // 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_x / 2, 0, cellsize_y });

            Kyokuritsuzu(src.Value, nrow - 1, ncol - 1, dst.Value);

            dst.Save();

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


        static void Kyokuritsuzu(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;

            for (int x = 1; x < ncol; x++)
            {
                for (int y = 1; y < nrow; y++)
                {
                    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 H22 = src[y, x];         if (H22 == -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 dx2 = Math.Pow(dx, 2);
                    double dy2 = Math.Pow(dy, 2);
                    double Kx = (H23 - H21) / (2 * dx);
                    double Ky = (H32 - H12) / (2 * dy);
                    double Kx2 = Math.Pow(Kx, 2);
                    double Ky2 = Math.Pow(Ky, 2);
                    double Kxx = (H21 + H23 - 2 * H22) / dx2;
                    double Kyy = (H12 + H32 - 2 * H22) / dy2;
                    double Kxy = (H31 - H13 - H33 + H11) / (2 * (dx2 + dy2));
                    double K = (Kxx * (1+ Ky2) + Kyy * (1+Kx2) -2 * Kx * Ky * Kxy) / (2 * Math.Pow(1+Kx2+Ky2, 3 / 2));
 
                    dst[y, x] = (K < 0 ? Math.Sqrt(Math.Abs(K)) * -1 : Math.Sqrt(K)) * 100;
                }
            }
            return;
        }
    }
}

準備

f:id:Chiakikun:20211230001444p:plain

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

 

 

 

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

DotSpatialを使って方位図も作成してみました。

概要

  • 前回に続いて、DEMから方位図を作成します。
  • 今回もラスタの読み書きにはDotSpatialを使用しました。
  • 傾斜量の逆正接関数によって方位を出しています。

f:id:Chiakikun:20211230000905p:plain

DEM

f:id:Chiakikun:20211230213840p:plain

方位図

f:id:Chiakikun:20211230214651p:plain

凡例

 

コード

using System;
using System.Runtime.InteropServices;

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

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

        static void Main(string[] args)
        {
            SetDllDirectory("DotSpatial.Data.Rasters.GdalExtensionのパス");

            string loadfilepath = "DEMファイル.tif";
            string savefilepath = "出力先.tif";

            // 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_x / 2, 0, cellsize_y });

            Shamenhoizu(src.Value, nrow - 1, ncol - 1, dst.Value);

            dst.Save();

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


        static void Shamenhoizu(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;

            for (int x = 1; x < ncol; x++)
            {
                for (int y = 1; y < nrow; y++)
                {
                    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 Hx = (H11 + H21 + H31 - (H13 + H23 + H33)) / (dx * 3);
                    double Hy = (H31 + H32 + H33 - (H11 + H12 + H13)) / (dy * 3);

                    if ((Hy == 0) && (Hx == 0))
                    {
                        dst[y, x] = -100;
                        continue;
                    }

                    double sita = Math.Atan(Hx / Hy) * 180.0 / Math.PI;
                    if (Hy < 0)
                    {
                        sita += 180;
                    }

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

16行目 DotSpatialをコンパイルしたときに作成されている「DotSpatial.Data.Rasters.GdalExtension」のフルパスを指定します。

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

準備

f:id:Chiakikun:20211230001444p:plain

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

 

 

 

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