C#とDotSpatialでGeoTiff読み書き
DotSpatialを使ってGeoTiffの読み書きをするプログラムのサンプルをご紹介いたします。
準備
①DotSpatialのソースをコンパイルします。詳しくは別のページで行っていますので興味のある方はそちらをご覧ください。
②ビルドが終わりました。
③Visual Studioを起動します。今回は上画像のプロジェクトを選択しました。
④プロジェクト名を設定後、このような画面になります。
⑤コンパイルしたDotSpatialをこのプログラムで使えるように、参照の追加をします。
今回は2つ追加しました。準備はここまでです。
書き出し
上の画像のようなGeoTiffを出力するプログラムのソースを紹介します。先ほど準備したプログラムに以下のソースを張り付けて実行すると、tifが作成されます。
ソース
using System.Runtime.InteropServices;
using DotSpatial.Data;
using DotSpatial.Data.Rasters.GdalExtension;
namespace ConsoleApp1
{
class Program
{
[DllImport("kernel32.dll", CharSet = CharSet.Auto)]
static extern bool SetDllDirectory(string lpPathName);
static void Main(string[] args)
{
// gdal等必要なdllをまとめて参照したいので
SetDllDirectory(@"D:\DotSpatial-master\Source\bin\Debug\Windows Extensions\DotSpatial.Data.Rasters.GdalExtension\gdal\x86");
string savefilepath = @"D:\test.tif";
int ncol = 3; // 水平方向ピクセル数
int nrow = 3; // 鉛直方向ピクセル数
double nodata = -9999;
double cellsize = 10;
double xllcorner = 20000;
double yllcorner = -40000;
// https://epsg.io/6677 から
string prj ="PROJCS[\"JGD2011 / Japan Plane Rectangular CS IX\","+
"GEOGCS[\"JGD2011\","+
"DATUM[\"Japanese_Geodetic_Datum_2011\"," +
"SPHEROID[\"GRS 1980\",6378137,298.257222101]," +
"TOWGS84[0, 0, 0, 0, 0, 0, 0]]," +
"PRIMEM[\"Greenwich\",0," +
"AUTHORITY[\"EPSG\",\"8901\"]]," +
"UNIT[\"degree\",0.0174532925199433," +
"AUTHORITY[\"EPSG\",\"9122\"]]," +
"AUTHORITY[\"EPSG\",\"6668\"]]," +
"PROJECTION[\"Transverse_Mercator\"]," +
"PARAMETER[\"latitude_of_origin\",36]," +
"PARAMETER[\"central_meridian\",139.833333333333]," +
"PARAMETER[\"scale_factor\",0.9999]," +
"PARAMETER[\"false_easting\",0]," +
"PARAMETER[\"false_northing\",0]," +
"UNIT[\"metre\",1," +
"AUTHORITY[\"EPSG\",\"9001\"]]," +
"AUTHORITY[\"EPSG\",\"6677\"]]";
// GeoTiff書き出し
GdalRasterProvider d = new GdalRasterProvider();
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[] { xllcorner, cellsize, 0, yllcorner, 0, -1 * cellsize });
dst.Value[0, 0] = 32; // H11
dst.Value[0, 1] = 64; // H12
dst.Value[0, 2] = 128; // H13
dst.Value[1, 0] = 16; // H21
dst.Value[1, 1] = 0; // H22
dst.Value[1, 2] = 1; // H23
dst.Value[2, 0] = 8; // H31
dst.Value[2, 1] = 4; // H32
dst.Value[2, 2] = 2; // H33
dst.Save();
return;
}
}
}
読み込み
先ほど書き出しのプログラムで作成したtifを読み込んで、ピクセルの値をコンソールに出力するプログラムのソースを紹介します。
ソース
using System;
using System.Runtime.InteropServices;
using DotSpatial.Data;
using DotSpatial.Data.Rasters.GdalExtension;
namespace ConsoleApp1
{
class Program
{
[DllImport("kernel32.dll", CharSet = CharSet.Auto)]
static extern bool SetDllDirectory(string lpPathName);
static void Main(string[] args)
{
// gdal等必要なdllをまとめて参照したいので
SetDllDirectory(@"D:\DotSpatial-master\Source\bin\Debug\Windows Extensions\DotSpatial.Data.Rasters.GdalExtension\gdal\x86");
string loadfilepath = @"D:\test.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];
Console.WriteLine("H11 = " + src.Value[0, 0]);
Console.WriteLine("H12 = " + src.Value[0, 1]);
Console.WriteLine("H13 = " + src.Value[0, 2]);
Console.WriteLine("H21 = " + src.Value[1, 0]);
Console.WriteLine("H22 = " + src.Value[1, 1]);
Console.WriteLine("H23 = " + src.Value[1, 2]);
Console.WriteLine("H31 = " + src.Value[2, 0]);
Console.WriteLine("H32 = " + src.Value[2, 1]);
Console.WriteLine("H33 = " + src.Value[2, 2]);
Console.ReadKey();
return;
}
}
}
最後までご覧頂き、ありがとうございました
PostGIS よく使うSQL文のメモ
- SQL文をファイルから実行する
- データベースの一覧表示
- PostGISのデータベースを作成する
- SQL文の実行時間を計測する
- 現在接続しているデータベース名を取得する
- データベースを削除する
- テーブル一覧
- テーブル作成
- テーブル名を変更する
- テーブルの行数を取得する
- テーブルのカラム名一覧
- テーブルを削除する
- テーブルを初期化する
- SRIDを変更する
- インデックスを作成する
- インデックスを削除する
- 点群のテキストを読み込む
このページでは自分がよく使うSQL文をメモしています。
SQL文をファイルから実行する
psql -U postgres -f ファイル名.sql
このコマンドを実行するためには環境変数の設定を済ませておく必要があります。
データベースの一覧表示
\l
PostGISのデータベースを作成する
CREATE DATABASE 〇〇; \c sample CREATE EXTENSION postgis;
1行目 〇〇というデータベースを作成しています。
2行目 〇〇に接続します。
3行目 現在のデータベースにPostGISの拡張をインストールしています。
SQL文の実行時間を計測する
\timing
このコマンドを実行した後、SQL文を実行すると、実行後に時間が出力されます。
現在接続しているデータベース名を取得する
SELECT current_database();
データベースを削除する
DROP DATABASE 〇〇;
テーブル一覧
\dt
テーブル作成
CREATE TABLE 〇〇
(
geom GEOMETRY(オブジェクトタイプ, SRID),
フィールド名 型
: :
: :
);
テーブル名を変更する
ALTER TABLE 〇〇 RENAME TO △△;
テーブルの行数を取得する
SELECT count(*) from 〇〇;
テーブルのカラム名一覧
\d 〇〇
テーブルを削除する
DROP TABLE 〇〇;
テーブルを初期化する
TRUNCATE 〇〇;
SRIDを変更する
SELECT UpdateGeometrySrID('〇〇', 'geom', SRID);
テーブル名とオブジェクトのフィールドには「'」を付けないとエラーになります。
インデックスを作成する
CREATE INDEX インデックス名 ON テーブル名 USING GIST (geom);
インデックスを削除する
DROP INDEX インデックス名;
点群のテキストを読み込む
スペースで区切られたx y 標高のテキストを読み込む場合...
CREATE TEMP TABLE tmp (
x double precision,
y double precision,
z double precision
);
\COPY tmp FROM 'points.txt' WITH CSV DELIMITER ' ';
CREATE TABLE pnt (
geom GEOMETRY(POINT, srid),
elev double precision
);
INSERT INTO
pnt (geom, elev)
SELECT
CONCAT ('SRID=srid;POINT(', x, ' ', y, ')') AS concated,
z
FROM
tmp;
1行目 tmpという一時的に使用するデータベースを作成しています。接続が切れると自動的に削除されます。
7行目 クライアント側に存在するファイル「points.txt」をテーブル「tmp」にコピーしています。テキストはスペース区切りです。
9行目 オブジェクトタイプがポイントのpntというテーブルを作成しています。
14行目 テーブル「tmp」のフィールド「x」の値xxxxと、フィールド「y」の値yyyyでオブジェクトを作成し、フィールド「z」の値と共に、「pnt」に投入しています。
PostGISのデータディレクトリを変更する
PostGISのインストール途中で「Failed to load SQL modules into the database cluster.」や「Spatial database creation failed.」とエラーが出る場合があり、データディレクトリを変更する必要がありましたので、変更手順を紹介します。
PostgreSQLのサービスを削除する
先に設定していたデータディレクトリを削除してから、PostgreSQLのインストーラーが設定したサービスを削除します。コントロールパネルから「システムとセキュリティ」→「管理ツール」→「サービス」と移動して「Postgresql-x64-10」が確認できます。
cmdを管理者権限で起動します。PowerShellではダメでした。
コマンドプロンプトで次のコマンドを実行します。
sc delete postgresql-x64-10
データディレクトリを作成する
データディレクトリを手動で作成します。ここはPowerShellでも構いません。initdbコマンドを使うのでPath設定はしておいてください。次のコマンドを実行します。
initdb -D C:\users\〇〇\PostgreSQL -W -U postgres -A password --no-locale -E UTF8
〇〇の部分は自身のユーザー名を入れてください。データベースのパスワードの入力を求められるのでパスワードを2回入力すると、データディレクトリが作成されます。
PostgreSQLのサービスを登録する
サービスを登録します。cmdを管理者権限で作成して、次のコマンドを実行します。
sc create postgresql-x64-10 binPath= "\"C:\Program Files\PostgreSQL\10\bin\pg_ctl.exe\" runservice -N \"postgresql-x64-10\" -D \"C:\Users\〇〇\PostgreSQL\"" start=auto password=************
