起伏量図を作成
概要
- 『幌延地域を対象とした10mグリッド数値標高モデルを用いた精密地形解析図の作成』を参考に、起伏量図を作成しました。
- ラスタの読み書きにはDotSpatialを使っています。
- ある範囲の高低差の大きさを表すことができます。
コード
using System;
using System.Runtime.InteropServices;
using DotSpatial.Data;
using DotSpatial.Data.Rasters.GdalExtension;
namespace Kifukuryozu
{
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:\DEM.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 });
Kifukuryo(src.Value, nrow, ncol, dst.Value);
dst.Save();
Console.WriteLine("終了しました。");
Console.ReadKey();
return;
}
static void Kifukuryo(IValueGrid src, int nrow, int ncol, IValueGrid dst)
{
int r = 5; // 窓サイズ
double[,] dsrc = new double[nrow, ncol];
for (int x = 0; x < ncol; x++)
for (int y = 0; y < nrow; y++)
dsrc[y, x] = src[y, x];
for (int x = 0; x < ncol; x++)
{
for (int y = 0; y < nrow; y++)
{
dst[y, x] = -9999;
// 窓領域内の最高標高、最低標高
double minelv = double.PositiveInfinity;
double maxelv = double.NegativeInfinity;
for (int m = x - r; m <= x + r; m++)
{
for (int n = y - r; n <= y + r; n++)
{
if ((dsrc[y, x] > -9999) && (m >= 0) && (n >= 0) && (m < ncol) && (n < nrow) && (dsrc[n, m] > -9999) )
{
minelv = Math.Min(minelv, dsrc[n, m]);
maxelv = Math.Max(maxelv, dsrc[n, m]);
}
}
}
if (double.IsInfinity(minelv) || double.IsInfinity(maxelv))
continue;
dst[y, x] = maxelv - minelv;
}
}
return;
}
}
}
準備
前回同様「DotSpatial.Data.dll」と「DotSpatial.Data.Rasters.GdalExtension.dll」を参照に追加してからソースをコンパイルしてください。
出力結果
起伏量が大きい場所から白、赤、黄、緑で色設定をしました。等高線間隔が狭い場所は白や赤になっているのがわかります。
最後までご覧いただき、ありがとうございました。
接峰面図・接谷面図を作成
概要
- 『幌延地域を対象とした10mグリッド数値標高モデルを用いた精密地形解析図の作成』を参考に、接峰面図・接谷面図を作成しました。
- ラスタの読み書きにはDotSpatialを使っています。
- プログラムは結構処理時間がかかります(ryzen5 2600で20分くらい)。もう少し何とかしたい...
コード
using System;
using System.Runtime.InteropServices;
using DotSpatial.Data;
using DotSpatial.Data.Rasters.GdalExtension;
namespace Seppomenzu
{
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:\DEM.tif";
string seppoupath = @"D:\接峰面.tif";
string sekkokupath = @"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 seppou = Raster.CreateRaster(seppoupath, null, ncol, nrow, 1, typeof(float), new[] { string.Empty });
seppou.NoDataValue = nodata;
seppou.ProjectionString = prj;
seppou.Bounds = new RasterBounds(nrow, ncol, new double[] { xllcenter - cellsize_x / 2, cellsize_x, 0, yllcenter - cellsize_y / 2, 0, cellsize_y });
IRaster sekkoku = Raster.CreateRaster(sekkokupath, null, ncol, nrow, 1, typeof(float), new[] { string.Empty });
sekkoku.NoDataValue = nodata;
sekkoku.ProjectionString = prj;
sekkoku.Bounds = new RasterBounds(nrow, ncol, new double[] { xllcenter - cellsize_x / 2, cellsize_x, 0, yllcenter - cellsize_y / 2, 0, cellsize_y });
Seppomen(src.Value, nrow, ncol, seppou.Value, sekkoku.Value);
seppou.Save();
sekkoku.Save();
Console.WriteLine("終了しました。");
Console.ReadKey();
return;
}
static void Seppomen(IValueGrid src, int nrow, int ncol, IValueGrid seppou, IValueGrid sekkoku)
{
int r = 5; // 窓サイズ
int cal = 10; // 計算回数
double[,] srcsep = new double[nrow, ncol];
double[,] srcsek = new double[nrow, ncol];
for (int x = 0; x < ncol; x++)
for (int y = 0; y < nrow; y++)
{
srcsep[y, x] = src[y, x];
srcsek[y, x] = src[y, x];
}
for (int i = 0; i < cal; i++)
{
double[,] dstsep = srcsep.Clone() as double[,];
double[,] dstsek = srcsek.Clone() as double[,];
for (int x = 0; x < ncol; x++)
{
for (int y = 0; y < nrow; y++)
{
dstsep[y, x] = -9999;
dstsek[y, x] = -9999;
if (srcsep[y, x] <= -9999)
continue;
// 窓領域内の平均
double cgridsep = 0;
double sumsep = 0;
double cgridsek = 0;
double sumsek = 0;
for (int m = x - r; m <= x + r; m++)
{
for (int n = y - r; n <= y + r; n++)
{
if ((m < 0) || (n < 0) || (m >= ncol) || (n >= nrow))
continue;
if ((srcsep[y, x] > -9999) && (srcsep[n, m] > -9999))
{
cgridsep += 1;
sumsep += srcsep[n, m];
}
if ((srcsek[y, x] > -9999) && (srcsek[n, m] > -9999))
{
cgridsek += 1;
sumsek += srcsek[n, m];
}
}
}
// 平均標高の方が元の標高値よりも低くなっているグリッドに元の標高値を代入
double avg = sumsep / cgridsep;
dstsep[y, x] = srcsep[y, x] <= avg ? avg : srcsep[y, x];
avg = sumsek / cgridsek;
dstsek[y, x] = srcsek[y, x] >= avg ? avg : srcsek[y, x];
}
}
srcsep = dstsep.Clone() as double[,];
srcsek = dstsek.Clone() as double[,];
}
for (int x = 0; x < ncol; x++)
for (int y = 0; y < nrow; y++)
{
seppou[y, x] = srcsep[y, x];
sekkoku[y, x] = srcsek[y, x];
}
}
}
}
準備
前回同様「DotSpatial.Data.dll」と「DotSpatial.Data.Rasters.GdalExtension.dll」を参照に追加してからソースをコンパイルしてください。
出力結果
下図の接峰面図では細かい谷が消え、尾根に接する部分ほど白く表示されています。接谷面は細かい尾根部分が削られて全体的に丸くなっています。カラーでは分かり辛いので黒白にしています。
最後までご覧いただき、ありがとうございました。
尾根谷度図を作成
概要
- 『幌延地域を対象とした10mグリッド数値標高モデルを用いた精密地形解析図の作成』を参考に、尾根谷度図を作成しました。
- ラスタの読み書きにはDotSpatialを使っています。
- 尾根谷度図の作成に必要な地上開度図と地下開度図はこちらで作成した図を使用しています。
コード
using System;
using System.Runtime.InteropServices;
using DotSpatial.Data;
using DotSpatial.Data.Rasters.GdalExtension;
namespace OneTanidozu
{
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 abovefilepah = @"D:\地上開度.tif";
string belowfilepah = @"D:\地下開度.tif";
string savefilepath = @"D:\保存先.tif";
SetDllDirectory(dllpath);
// GeoTiff読み込み
GdalRasterProvider d = new GdalRasterProvider();
IRaster above = d.Open(abovefilepah);
IRaster below = d.Open(belowfilepah);
int ncol = above.NumColumns;
int nrow = above.NumRows;
int band_num = above.NumBands;
string prj = above.ProjectionString;
double nodata = above.NoDataValue;
double[] pGT = above.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 });
OneTanidozu(above.Value, below.Value, nrow, ncol, dst.Value);
dst.Save();
Console.WriteLine("終了しました。");
Console.ReadKey();
return;
}
static void OneTanidozu(IValueGrid above, IValueGrid below, int nrow, int ncol, IValueGrid dst)
{
for (int x = 0; x < ncol; x++)
{
for (int y = 0; y < nrow; y++)
{
dst[y, x] = -9999;
double az = above[y, x];
double bz = below[y, x];
if (az <= -9999 || bz <= -9999)
continue;
dst[y, x] = (az - bz) / 2.0;
}
}
return;
}
}
}
準備
前回同様「DotSpatial.Data.dll」と「DotSpatial.Data.Rasters.GdalExtension.dll」を参照に追加してからソースをコンパイルしてください。
出力結果
下図は起伏が緩い場所の尾根谷度図です。谷が黒く、尾根は白く表示されるようにしています。地形の変わり目がくっきり見やすくなったと思います。
最後までご覧いただき、ありがとうございました。
