水系密度図を作成
概要
- 『幌延地域を対象とした10mグリッド数値標高モデルを用いた精密地形解析図の作成』を参考に、前回作成した尾根・谷図を使って水系密度図を作成してみました。
- 今回もラスタの読み書きにはDotSpatialを使用しています。
- 斜面の凹凸具合を見ることができます。
コード
using System; using System.Runtime.InteropServices; using DotSpatial.Data; using DotSpatial.Data.Rasters.GdalExtension; namespace Suikeimitsudozu { class Program { [DllImport("kernel32.dll", CharSet = CharSet.Auto)] static extern bool SetDllDirectory(string lpPathName); static void Main(string[] args) { string dllpath = @"D:\DotSpatial-master\Source\bin\Debug\Windows Extensions\DotSpatial.Data.Rasters.GdalExtension\gdal\x86"; string loadfilepath = @"D:\尾根谷図.tif"; // 尾根谷図の出力を使う string savefilepath = @"D:\保存先.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 }); Suikeimitsudo(src.Value, nrow, ncol, dst.Value); dst.Save(); Console.WriteLine("終了しました。"); Console.ReadKey(); return; } static bool isTani(double value) { return value >= 1 ? true : false; } static void Suikeimitsudo(IValueGrid gsrc, int nrow, int ncol, IValueGrid dst) { int r = 10; double[,] src = new double[nrow, ncol]; for (int x = 0; x < ncol; x++) for (int y = 0; y < nrow; y++) { if (src[y, x] <= -9999) src[y, x] = 0; else src[y, x] = gsrc[y, x]; } for (int x = 0; x < ncol; x++) { for (int y = 0; y < nrow; y++) { dst[y, x] = -9999; int count = 0; int fluct_x = 0, fluct_y = 0; for (int x_axis = 0; x_axis <= r; x_axis++) { int y_axis = (int)Math.Truncate(Math.Sqrt(Math.Pow(r, 2) - Math.Pow(x_axis, 2))); for (int i = 0; i <= y_axis; i++) { // 第1象限 fluct_x = x_axis; fluct_y = i; if ((x + fluct_x >= 0) && (y + fluct_y >= 0) && (x + fluct_x < ncol) && (y + fluct_y < nrow)) { if (isTani(src[y + fluct_y, x + fluct_x])) count++; } // 第2象限 fluct_x = x_axis * (-1); fluct_y = i; if(fluct_x != 0) // 第1象限で計算済だから { if ((x + fluct_x >= 0) && (y + fluct_y >= 0) && (x + fluct_x < ncol) && (y + fluct_y < nrow)) { if (isTani(src[y + fluct_y, x + fluct_x])) count++; } } // 第3象限 fluct_x = x_axis * (-1); fluct_y = i * (-1); if (fluct_y != 0) // 第2象限で計算済だから { if ((x + fluct_x >= 0) && (y + fluct_y >= 0) && (x + fluct_x < ncol) && (y + fluct_y < nrow)) { if (isTani(src[y + fluct_y, x + fluct_x])) count++; } } // 第4象限 fluct_x = x_axis; fluct_y = i * (-1); if ((fluct_y != 0) && (fluct_x != 0)) // 第1、第3象限で計算済だから { if ((x + fluct_x >= 0) && (y + fluct_y >= 0) && (x + fluct_x < ncol) && (y + fluct_y < nrow)) { if (isTani(src[y + fluct_y, x + fluct_x])) count++; } } } } dst[y, x] = count; } } return; } } }
準備
前回同様「DotSpatial.Data.dll」と「DotSpatial.Data.Rasters.GdalExtension.dll」を参照に追加してからソースをコンパイルしてください。
出力結果
下図の左に谷が多く存在するので、その部分が水系密度図では赤くなっています。
最後までご覧いただき、ありがとうございました。