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」を参照に追加してからソースをコンパイルしてください。
最後までご覧いただき、ありがとうございました。