方位図を作成する
概要
- 『幌延地域を対象とした10mグリッド数値標高モデルを用いた精密地形解析図の作成』を参考に、方位図を作成しました。
- 今回もラスタの読み書きには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) { 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 }); Shamenhoizu(src.Value, nrow, ncol, 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; 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 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; } } }
58、59行目の値は、参考資料のP3に載っていた縦方向、横方向の1グリッド分の距離です。
準備
前回同様「DotSpatial.Data.dll」と「DotSpatial.Data.Rasters.GdalExtension.dll」を参照に追加してからソースをコンパイルしてください。
凡例について
上の画像で使っている凡例はこちらに置いてあります。よろしければお使いください。
出力結果
円錐形の山だと方位の色がきれいに表示されるだろうと思い、長野県の蓼科山と鹿児島県の開聞岳の方位図を出力してみました。
最後までご覧いただき、ありがとうございました。