[QGIS & Gdal]QGISにてGeotiffファイルのピクセルサイズが一定でない? ~その2~

2019/7/16 C#, DotSpatial, gdal, GIS, QGIS, プログラミング

前回の記事で「QGISでGeotiffファイルを拡大するとピクセルサイズが一定でない」と書きましたが、その発生条件について追試してみたのでその結果を。

追試内容は

  1. 平面直角座標系(x,y座標/m単位)のメッシュデータだとどうなるのか？・・・メッシュサイズ：250mの場合
  2. 上記の場合でメッシュサイズに端数が出る場合だとどうなるのか？・・・メッシュサイズ：800/3＝266.666…mの場合

の二ケースです。

早速、試してみました。
座標系はUTM54帯(JGD2000、EPSG:3100)としました。
まずは250mメッシュの場合。

250mメッシュの場合
点はCSVファイルに記録したx,y座標によるもの(ラスタの座標と同じ)
ラスタの値はランダム値に基づくグレースケール


見ての通り、ずれていません。

次に266.666…mメッシュの場合。

266.666…mメッシュの場合

こちらもずれていません。
どうやら、緯度・経度の座標系で、メッシュサイズの数値が小さい(緯度経度の250mメッシュの場合、7.5秒＝0.00208333…度)場合に何らかの誤差が発生しているようです。

QGISをご利用の方で、緯度経度単位で作成されたラスタデータを使われる方はご注意ください。
当方の環境：QGIS2.18.11、64bit、Windows10

最後に、266mメッシュのGeotiffを作った時のコードを。
下記のメソッドは省略しましたので、それぞれのリンク先をご参照ください。
InitGDAL():GDALを使用するためのメソッド
GetRandom():一様乱数を発生させるメソッド

const int iW = 40;
const int iH = 40;
const string sOutCsv = @"C:\…\Mesh_xy266.csv";
const string sOutTif = @"C:\…\Mesh_xy266.tif";
const double dLon0 = 318000d;
const double dLat0 = 3925000d;
const double dX = 800d / 3; //1pxあたりの幅　：800/3＝266.66…
const double dY = 800d / 3; //1pxあたりの高さ：800/3＝266.66…
const string sPrj = "PROJCS[\"JGD2000 / UTM zone 54N\",GEOGCS[\"JGD2000\",DATUM[\"Japanese_Geodetic_Datum_2000\",SPHEROID[\"GRS 1980\",6378137,298.257222101,AUTHORITY[\"EPSG\",\"7019\"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY[\"EPSG\",\"6612\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4612\"]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",141],PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",500000],PARAMETER[\"false_northing\",0],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH],AUTHORITY[\"EPSG\",\"3100\"]]";

static void Main(string[] args)
{
    Encoding sjis = Encoding.GetEncoding(932);
    InitGdal(); //GDALを初期化する(環境変数の設定など)

    double[] dV = new double[iH * iW];
    using(StreamWriter sw = new StreamWriter(sOutCsv, false, Encoding.GetEncoding(932)))
    {
        int i = 0;
        sw.WriteLine("Lon,Lat,Value");
        for (int iY = 0; iY < iH; iY++)
        {
            for (int iX = 0; iX < iW; iX++)
            {
                dV[i] = GetRandom();
                sw.WriteLine("{0},{1},{2}", dLon0 + dX / 2 + dX * iX, dLat0 + dY / 2 + dY * iY, dV[i]);
                i++;
            }
        }
    }

    GDAL.Driver drv = GDAL.Gdal.GetDriverByName("GTiff");
    
    //出力ファイル：ファイル名・幅・高さ・バンド数・ピクセルのデータ型、オプション(null)
    using (GDAL.Dataset dsOut = drv.Create(sOutTif, iW, iH, 1, GDAL.DataType.GDT_Float64, null))
    {
        //位置情報(Xmin, ΔX/px, 0, Ymax, 0, ΔY/px)
        double[] dGeo = new double[] { dLon0, dX, 0d, dLat0, 0, dY };
        dsOut.SetGeoTransform(dGeo);

        using (GDAL.Band bd = dsOut.GetRasterBand(1))
        {
            //データ書き込み
            bd.WriteRaster(0, 0, iW, iH, dV, iW, iH, 0, 0);
            bd.FlushCache();
        }
        dsOut.SetProjection(sPrj);
        dsOut.FlushCache();
    }
}

//以下省略(InitGidal、GetRandom)