DotSpatialでGeoTiff読み書き

DotSpatialを使ってGeoTiffの読み書きを行うサンプルを残しておきます。

DotSpatial版の他に、DotSpatialをインストールすると付いてくるgdal_csharp.dllを使った場合のソースも置いてます。


DotSpatial版

 [DllImport("kernel32.dll", CharSet = CharSet.Auto)]
  static extern bool SetDllDirectory(string lpPathName);

  static void Main(string[] args)
  {
    SetDllDirectory(gdal_wrap.dllが置いてあるフォルダパス。例えば@"~\DotSpatial-master\Source\bin\x64\Release\gdal\x64"とか);

    string loadfilepath = 読み込むファイルパス;
    string savefilepath = 書き込みファイルパス;

    // 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];
    // pGT[2] # 回転 今回は使わない
    double yllcenter = pGT[3];
    // pGT[4] # 回転 今回は使わない
    double cellsize_y = pGT[5];

    IRaster dst = Raster.CreateRaster(savefilepath, src.DriverCode, ncol - 2, nrow - 2, 1, src.DataType, new[] { string.Empty });
    dst.NoDataValue = nodata;
    dst.ProjectionString = prj;
    dst.Bounds = new RasterBounds(nrow - 2, ncol - 2,
        new double[] { xllcenter + cellsize_x / 2, cellsize_x, 0, yllcenter + cellsize_y / 2, 0, cellsize_y });

    Hoi(src.Value, nodata, xllcenter, yllcenter, cellsize_x, nrow, ncol, dst.Value);

    dst.Save();
  }

1,2,6行目 「gdal_wrap.dll」を読み込むためにこの記述が必要です。『参照の追加』ではこのdllを読み込めないのです。

12,13行目 GeoTiffを読み込んでいます。

14~26行目 読み込んだGeoTiffの情報を取得しています。

28~32行目 書き込み用GeoTiffを作成しています。このサンプルでは、読み込んだGeoTiffより一回り小さいGeoTiffを作成するので、サイズを基画像の縦横-2(31行目)、原点を1ピクセル分内側(32行目)に設定しました。他にGeoTiffの原点、ピクセルのサイズを設定しています。

34行目 読み込んだGeoTiffを基に、斜面方位を計算して、書き込み用GeoTiffのピクセルデータを作成しています。(今回は説明を省きます)

36行目 書き込み用GeoTiffを保存しています。

 

ソースはこちらです。ソース見れない状態でした。すいませんでした。


gdal_csharp版

  [DllImport("kernel32.dll", CharSet = CharSet.Auto)]
  static extern bool SetDllDirectory(string lpPathName);

  static void Main(string[] args)
  {
    SetDllDirectory(gdal_wrap.dllが置いてあるフォルダパス。例えば@"~\DotSpatial-master\Source\bin\x64\Release\gdal\x64"とか);

    string loadfilepath = 読み込むファイルパス;
    string savefilepath = 書き込みファイルパス;

    Gdal.AllRegister();

    // GeoTiff読み込み
    Dataset ds = Gdal.Open(loadfilepath, Access.GA_Update);
    int ncol = ds.RasterXSize; // 水平方向ピクセル数
    int nrow = ds.RasterYSize; // 鉛直方向ピクセル数
    int band_num = ds.RasterCount; // バンド数
    string prj = ds.GetProjection();

    double nodata;
    int hasval;
    ds.GetRasterBand(1).GetNoDataValue(out nodata, out hasval);

    double[] pGT = new double[6];
    ds.GetGeoTransform(pGT);
    double xllcorner = pGT[0];
    double cellsize_x = pGT[1];
    // pGT[2] # 回転 今回は使わない
    double yllcorner = pGT[3];
    // pGT[4] # 回転 今回は使わない
    double cellsize_y = pGT[5];

    double[][] dVals = new double[band_num][];
    for (int i = 1; i <= band_num; i++)
    {
        dVals[i - 1] = new double[ncol * nrow];
        ds.GetRasterBand(i).ReadRaster(0, 0, ncol, nrow, dVals[i - 1], ncol, nrow, 0, 0);
    }

    double xllcenter = xllcorner + cellsize_x / 2;
    double yllcenter = yllcorner + cellsize_y / 2;
    double[] output = Hoi(dVals[0], nodata, xllcenter, yllcenter, cellsize_x, nrow, ncol);

    // GeoTiff書き出し
    ds = null;
    Driver imgDriver = Gdal.GetDriverByName("GTiff");
    ds = imgDriver.Create(savefilepath, ncol - 2, nrow - 2, 1, DataType.GDT_Float32, new string[] { });
    ds.SetProjection(prj);
    ds.SetGeoTransform(new double[] { xllcorner + cellsize_x, cellsize_x, 0, yllcorner - cellsize_x, 0, -1 * cellsize_x });
    Band rband = ds.GetRasterBand(1);
    rband.WriteRaster(0, 0, ncol - 2, nrow - 2, output, ncol - 2, nrow - 2, 0, 0);
    rband.SetNoDataValue(nodata);
    ds.FlushCache();

    return;
  }

11行目 GDALのドライバを登録しています。アプリケーションの開始時に1回呼び出す必要があります。

14行目 GeoTiffを読み込んでいます。

15~38行目 読み込んだTiffの情報を取得しています。DotSpatialと違い、原点は画像の端っこの座標を取得します(26行目と29行目)。33行目で画像の各ピクセルデータをdoubleの配列に格納しています。

40~42行目 原点ピクセルの中心点を取得して、42行目で斜面方位を計算してdoubleの配列を取得しています。

46,47行目 書き込む空のGeoTiffを作成しています。サイズは今回も基画像より一回り小さいサイズです。

49行目 書き込むGeoTiffの原点やピクセルサイズを設定しています。

51,53行目 42行目で作成したピクセルデータをディスクに書き込んでいます。(52行目は上に移した方が見栄え良いと思います。)

 

ソースはこちらです。ソース見れない状態でした。すいませんでした。

最後までご覧頂き、ありがとうございました

PostGISメモ

このページでは自分がよく使うSQL文をメモしています。


SQL文をファイルから実行する

psql -U postgres -f ファイル名.sql

このコマンドを実行するためには環境変数の設定を済ませておく必要があります。


PostGISのデータベースを作成する

CREATE DATABASE sample;
\c sample
CREATE EXTENSION postgis;
\l

1行目 sampleというデータベースを作成しています。

2行目 sampleに接続します。

3行目 現在のデータベースにPostGISの拡張をインストールしています。

4行目 データベースの一覧を表示しています。


点群のテキストを読み込む

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, 4326),
elev double precision
);

INSERT INTO pnt (geom, elev) SELECT CONCAT ('SRID=4326;POINT(', x, ' ', y, ')') AS concated, z FROM tmp;

1行目 tmpという一時的に使用するデータベースを作成しています。接続が切れると自動的に削除されます。

8行目 クライアント側に存在するファイル「points.txt」をテーブル「tmp」にコピーしています。テキストはスペース区切りです。

10行目 オブジェクトタイプがポイントのpntというテーブルを作成しています。

15行目 テーブル「tmp」のフィールド「x」の値xxxxと、フィールド「y」の値yyyyを使って、文字列「SRID=4326;POINT(xxxx yyyy)」を作成し、concatedに入れています。

テーブル「pnt」のフィールド「geom」にconcatedの値を、フィールド「elev」にテーブル「tmp」のフィールド「z」の値を入れています。


ポリゴンに重なるポイントを取得する

CREATE TABLE intpnt (LIKE pnt);
INSERT INTO intpnt (geom, elev) (SELECT geom, elev  FROM pnt WHERE ST_Intersects(geom, (SELECT geom FROM mesh where コード='53381465')));

1行目 構造はpntと同じintpntというデータベースを作成しています。

2行目 テーブル「intpnt」に、テーブル「mesh」のフィールド「コード」の値が53381465のレコードのオブジェクトに重なるテーブル「pnt」のオブジェクトを追加します。


テーブル名を変更する

ALTER TABLE 〇〇〇 RENAME TO △△△;

間違えてテーブル名の先頭に数字を使ってしまったとき、〇〇〇はダブルクォーテーションで囲まないと、syntax errorになります。

PostGISのデータディレクトリを変更する


PostGISのインストール途中で「Failed to load SQL modules into the database cluster.」や「Spatial database creation failed.」とエラーが出る場合があり、データディレクトリを変更する必要がありましたので、変更手順を紹介します。


PostgreSQLのサービスを削除する

f:id:Chiakikun:20190815233634p:plain

先に設定していたデータディレクトリを削除してから、PostgreSQLインストーラーが設定したサービスを削除します。コントロールパネルから「システムとセキュリティ」→「管理ツール」→「サービス」と移動して「Postgresql-x64-10」が確認できます。

 

f:id:Chiakikun:20190815234039p:plain

cmdを管理者権限で起動します。PowerShellではダメでした。

コマンドプロンプトで次のコマンドを実行します。

sc delete postgresql-x64-10

 


データディレクトリを作成する

f:id:Chiakikun:20190815234342p:plain

データディレクトリを手動で作成します。ここはPowerShellでも構いません。initdbコマンドを使うのでPath設定はしておいてください。次のコマンドを実行します。

initdb -D C:\users\〇〇\PostgreSQL -W -U postgres -A password --no-locale -E UTF8

〇〇の部分は自身のユーザー名を入れてください。データベースのパスワードの入力を求められるのでパスワードを2回入力すると、データディレクトリが作成されます。


PostgreSQLのサービスを登録する

f:id:Chiakikun:20190815235038p:plain

サービスを登録します。cmdを管理者権限で作成して、次のコマンドを実行します。

sc create postgresql-x64-10 binPath= "\"C:\Program Files\PostgreSQL\10\bin\pg_ctl.exe\" runservice -N \"postgresql-x64-10\" -D \"C:\Users\me\PostgreSQL\"" start=auto password=************