DotSpatialを使って遷急線・遷緩線図と尾根・谷線図も作成してみました。
概要
コード
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; } } }
準備
前回同様「DotSpatial.Data.dll」と「DotSpatial.Data.Rasters.GdalExtension.dll」を参照に追加してからソースをコンパイルしてください。
最後までご覧いただき、ありがとうございました。
DotSpatialを使って曲率図も作成してみました。
概要
- 前回に続いて、DEMから曲率図を作成します。
- 今回もラスタの読み書きにはDotSpatialを使用しました。
- 今月作成したページは、画像とコードを入れ替えてるだけの手抜き仕様。
コード
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; } } }
準備
前回同様「DotSpatial.Data.dll」と「DotSpatial.Data.Rasters.GdalExtension.dll」を参照に追加してからソースをコンパイルしてください。
最後までご覧いただき、ありがとうございました。
DotSpatialを使って方位図も作成してみました。
概要
- 前回に続いて、DEMから方位図を作成します。
- 今回もラスタの読み書きにはDotSpatialを使用しました。
- 傾斜量の逆正接関数によって方位を出しています。
コード
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グリッド分の距離です。
準備
前回同様「DotSpatial.Data.dll」と「DotSpatial.Data.Rasters.GdalExtension.dll」を参照に追加してからソースをコンパイルしてください。
最後までご覧いただき、ありがとうございました。