ラプラシアン図とエッジ図を作成してみました。
概要
- 『幌延地域を対象とした10mグリッド数値標高モデルを用いた精密地形解析図の作成』を参考に、ラプラシアン図とエッジ図を作成しました。
- 今回もラスタの読み書きにはDotSpatialを使用しました。
コード
using System;
using System.Runtime.InteropServices;
using DotSpatial.Data;
using DotSpatial.Data.Rasters.GdalExtension;
namespace Laplacianzu
{
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 savelappath = @"D:\ラプラシアン.tif";
string saveedgpath = @"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 lap = Raster.CreateRaster(savelappath, null, ncol, nrow, 1, typeof(float), new[] { string.Empty });
lap.NoDataValue = nodata;
lap.ProjectionString = prj;
lap.Bounds = new RasterBounds(nrow, ncol, new double[] { xllcenter - cellsize_x / 2, cellsize_x, 0, yllcenter - cellsize_y / 2, 0, cellsize_y });
IRaster edg = Raster.CreateRaster(saveedgpath, null, ncol, nrow, 1, typeof(float), new[] { string.Empty });
edg.NoDataValue = nodata;
edg.ProjectionString = prj;
edg.Bounds = new RasterBounds(nrow, ncol, new double[] { xllcenter - cellsize_x / 2, cellsize_x, 0, yllcenter - cellsize_y / 2, 0, cellsize_y });
Laplacianzu(src.Value, nrow, ncol, lap.Value, edg.Value);
lap.Save();
edg.Save();
Console.WriteLine("終了しました。");
Console.ReadKey();
return;
}
static void Laplacianzu(IValueGrid src, int nrow, int ncol, IValueGrid lap, IValueGrid edg)
{
double dx = 8.763;
double dy = 12.348;
for (int x = 0; x < ncol; x++)
{
for (int y = 0; y < nrow; y++)
{
lap[y, x] = -9999;
edg[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 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 I = ((H21 + H23) - 2 * H22) / Math.Pow(dx, 2) + ((H12 + H32) - 2 * H22) / Math.Pow(dy, 2);
lap[y, x] = I;
edg[y, x] = Math.Abs(I);
}
}
return;
}
}
}
準備
前回同様「DotSpatial.Data.dll」と「DotSpatial.Data.Rasters.GdalExtension.dll」を参照に追加してからソースをコンパイルしてください。
出力結果
下図の左下は地形の変化が少ないのでラプラシアン図、エッジ図も白い部分が多く、右上は変化が多いので色が出ている部分が多くなっています。
『幌延地域を対象とした10mグリッド数値標高モデルを用いた精密地形解析図の作成』を参考に各種フィルタ作成をしてきましたが、これで終わりです。ご覧いただきありがとうございました。
高度分散量異常量図を作成
概要
- 『幌延地域を対象とした10mグリッド数値標高モデルを用いた精密地形解析図の作成』を参考に、高度分散量異常量図を作成しました。
- ラスタの読み書きにはDotSpatialを使っています。
- ここで使用している高度分散量図は以前作成したものを使用しています。
コード
using System;
using System.Runtime.InteropServices;
using System.Linq;
using DotSpatial.Data;
using DotSpatial.Data.Rasters.GdalExtension;
using System.Collections.Generic;
namespace KodoBunsanryoIjouzu
{
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 loaddempath = @"D:\DEM.tif";
string loadkodopath = @"D:\高度分散量図.tif";
string savefilepath = @"D:\保存先.tif";
SetDllDirectory(dllpath);
// GeoTiff読み込み
GdalRasterProvider d = new GdalRasterProvider();
IRaster dem = d.Open(loaddempath);
IRaster kod = d.Open(loadkodopath);
int ncol = dem.NumColumns;
int nrow = dem.NumRows;
int band_num = dem.NumBands;
string prj = dem.ProjectionString;
double nodata = dem.NoDataValue;
double[,] avgKodoBunsan = AvgKodoBunsanryo(dem.Value, kod.Value, nrow, ncol);
double [] keisu = GetKinjiKeisu(avgKodoBunsan);
double[] pGT = dem.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 });
Ijouzu(dem.Value, kod.Value, nrow, ncol, dst.Value, keisu);
dst.Save();
Console.WriteLine("終了しました。");
Console.ReadKey();
return;
}
static double[,] AvgKodoBunsanryo(IValueGrid dem, IValueGrid kodo, int nrow, int ncol)
{
List[] avgs = new List[29]; // 幌延町だとこのくらい必要
for (int i = 0; i < avgs.Length; i++)
avgs[i] = new List();
for (int x = 0; x < ncol; x++)
{
for (int y = 0; y < nrow; y++)
{
if (dem[y, x] <= -9999)
continue;
int idx = (int)dem[y, x] / 25;
avgs[idx].Add(kodo[y, x]);
}
}
double[,] ret = new double[avgs.Length, 2];
for (int i = 0; i < avgs.Length; i++)
{
ret[i, 0] = i * 25;
ret[i, 1] = avgs[i].Sum() / avgs[i].Count;
}
return ret;
}
static public double[] GetKinjiKeisu(double[,] pnts)
{
double[,] Det_a = new double[4, 4];
double[,] Inv_a = new double[4, 4];
double[,] Det_b = new double[4, 1];
double[,] Det_c = new double[4, 1];
double buf;
int m = 4;
for (int i = 1; i < pnts.GetLength(0) - 2; i++) // 資料の見た目に合わせるために調整した
{
double x6 = Math.Pow(pnts[i, 0], 6);
double x5 = Math.Pow(pnts[i, 0], 5);
double x4 = Math.Pow(pnts[i, 0], 4);
double x3 = Math.Pow(pnts[i, 0], 3);
double x2 = Math.Pow(pnts[i, 0], 2);
double x = pnts[i, 0];
Det_a[0, 0] += x6; Det_a[0, 1] += x5; Det_a[0, 2] += x4; Det_a[0, 3] += x3;
Det_a[1, 0] += x5; Det_a[1, 1] += x4; Det_a[1, 2] += x3; Det_a[1, 3] += x2;
Det_a[2, 0] += x4; Det_a[2, 1] += x3; Det_a[2, 2] += x2; Det_a[2, 3] += x;
Det_a[3, 0] += x3; Det_a[3, 1] += x2; Det_a[3, 2] += x; Det_a[3, 3] = 26;
Det_b[0, 0] += x3 * pnts[i, 1];
Det_b[1, 0] += x2 * pnts[i, 1];
Det_b[2, 0] += x * pnts[i, 1];
Det_b[3, 0] += pnts[i, 1];
}
for (int i = 0; i < m; i++)
for (int j = 0; j < m; j++)
Inv_a[i, j] = (i == j) ? 1.0 : 0.0;
for (int i = 0; i < m; i++)
{
buf = 1 / Det_a[i, i];
for (int j = 0; j < m; j++)
{
Det_a[i, j] *= buf;
Inv_a[i, j] *= buf;
}
for (int j = 0; j < m; j++)
if (i != j)
{
buf = Det_a[j, i];
for (int k = 0; k < m; k++)
{
Det_a[j, k] -= Det_a[i, k] * buf;
Inv_a[j, k] -= Inv_a[i, k] * buf;
}
}
}
for (int i = 0; i < 4; i++)
for (int j = 0; j < 1; j++)
for (int k = 0; k < 4; k++)
Det_c[i, j] += Inv_a[i, k] * Det_b[k, j];
double[] result = new double[4];
result[0] = Det_c[0, 0];
result[1] = Det_c[1, 0];
result[2] = Det_c[2, 0];
result[3] = Det_c[3, 0];
return result;
}
static void Ijouzu(IValueGrid dem, IValueGrid kodo, int nrow, int ncol, IValueGrid dst, double []keisu)
{
double a = keisu[0];
double b = keisu[1];
double c = keisu[2];
double d = keisu[3];
for (int x = 0; x < ncol; x++)
{
for (int y = 0; y < nrow; y++)
{
if (dem[y, x] <= -9999) continue;
double kinji = a * Math.Pow(dem[y, x], 3) + b * Math.Pow(dem[y, x], 2) + c * dem[y, x] + d;
dst[y, x] = kodo[y, x] - kinji;
}
}
return;
}
}
}
メモ コードの37行目で平均高度分散量（下図の青色線）、38行目で高度分散量近似曲線（下図の橙色線）を取得しています。参考資料に合わせたかったのでプログラムを調整しています。
準備
前回同様「DotSpatial.Data.dll」と「DotSpatial.Data.Rasters.GdalExtension.dll」を参照に追加してからソースをコンパイルしてください。
出力結果
近似曲線から特に大きく外れている箇所の高度分散量図と高度分散量異常量図を並べてみました。何で大きく外れることとなったのかは...わかりません。
最後までご覧いただき、ありがとうございました。
高度分散量図を作成
概要
- 『幌延地域を対象とした10mグリッド数値標高モデルを用いた精密地形解析図の作成』を参考に、高度分散量図を作成しました。
- ラスタの読み書きにはDotSpatialを使っています。
- 起伏量図よりも地形がわかり易く表されています。
コード
using System;
using System.Runtime.InteropServices;
using DotSpatial.Data;
using DotSpatial.Data.Rasters.GdalExtension;
using System.Collections.Generic;
namespace KodoBunsanryozu
{
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 });
KodoBunsanryo(src.Value, nrow, ncol, dst.Value);
dst.Save();
Console.WriteLine("終了しました。");
Console.ReadKey();
return;
}
static void KodoBunsanryo(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;
if (dsrc[y, x] <= -9999)
continue;
// 窓領域内の平均
double cgrid = 0;
double sum = 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) && (dsrc[n, m] > -9999))
{
sum += dsrc[n, m];
cgrid += 1;
}
}
}
double avg = sum / cgrid;
// 偏差
double dev = 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) && (dsrc[n, m] > -9999))
dev += Math.Pow(dsrc[n, m] - avg, 2);
}
dev = Math.Sqrt(dev / cgrid);
dst[y, x] = dev;
}
}
return;
}
}
}
準備
前回同様「DotSpatial.Data.dll」と「DotSpatial.Data.Rasters.GdalExtension.dll」を参照に追加してからソースをコンパイルしてください。
出力結果
起伏量図だと地形がギザギザ？した感じに見えますが、高度分散量図だと地形がわかり易くなりました。色設定は起伏量図と同じように設定しています。
最後までご覧いただき、ありがとうございました。
