前回はGDAL/OGRを使って点(Point)を作図するという記事を書きましたが、今回はポリゴンを作図してみたいと思います。
基本的な流れは前回と同じなのですが、Polygonのジオメトリ(wkbPolygon)はポリゴンを構成するリング(wkbLinearRing)から構成されているのがポイントの場合とは異なるところです。
また、穴あきポリゴンも、二つのリングをPolygonのジオメトリに追加してあげることにより作図することが出来ます。
サンプルソースでは、50~55行目で外側リング、58~63行目で内側リングを作成し、これらのリングを元に66~68行目でポリゴンを作図しています。
//これが書いてある前提
using System.IO;
using OGR = OSGeo.OGR;
using OSR = OSGeo.OSR;
/*中略*/
const string sGdalDir = @"C:\Program Files\GDAL"; //GDALのインストールディレクトリ
const string sShp = @"D:\hogehoge\2_poly.shp";
//投影法:OGC WKTフォーマット
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)
{
#region "OGRの準備(環境変数Pathの設定→Ogr.RegisterAll())"
string sVal = Environment.GetEnvironmentVariable("PATH") + ";" +
sGdalDir + ";" +
Path.Combine(sGdalDir, "csharp") + ";" +
Path.Combine(sGdalDir, "gdal-data") + ";" +
Path.Combine(sGdalDir, "gdalplugins") + ";" +
Path.Combine(sGdalDir, "projlib");
Environment.SetEnvironmentVariable("PATH", sVal);
OGR.Ogr.RegisterAll();
#endregion
//Shpを作成する
//Driver → DataSource → Layer
OGR.Driver drv = OGR.Ogr.GetDriverByName("ESRI Shapefile");
using (OGR.DataSource ds = drv.CreateDataSource(sShp, null))
using (OGR.Layer layer = ds.CreateLayer("layer", new OSR.SpatialReference(sPrj), OGR.wkbGeometryType.wkbPolygon, null))
{
OGR.FeatureDefn fetDef = layer.GetLayerDefn(); //レイヤと地物を結びつける?
//属性テーブルにフィールド追加
OGR.FieldDefn fdStr = new OGR.FieldDefn("fStr", OGR.FieldType.OFTString); //文字列型
layer.CreateField(fdStr, 1);
OGR.FieldDefn fldDef = new OGR.FieldDefn("fX", OGR.FieldType.OFTReal); //実数型
layer.CreateField(fldDef, 1);
for (long l = 0; l < 100; l++) //地物を100個作成
{
double dX = 10 * (l % 10); double dY = 10 * (l / 10); //基準となる座標
//ポリゴン作図 参考→https://pcjericks.github.io/py-gdalogr-cookbook/geometry.html
//外側の線
OGR.Geometry geoOutR = new OGR.Geometry(OGR.wkbGeometryType.wkbLinearRing);
geoOutR.AddPoint(dX, dY, 0d);
geoOutR.AddPoint(dX + 5d, dY, 0d);
geoOutR.AddPoint(dX + 5d, dY + 5d, 0d);
geoOutR.AddPoint(dX, dY + 5d, 0d);
geoOutR.CloseRings(); //閉じる
//内側の線
OGR.Geometry geoInR = new OGR.Geometry(OGR.wkbGeometryType.wkbLinearRing);
geoInR.AddPoint(dX + 1, dY + 1, 0d);
geoInR.AddPoint(dX + 3, dY + 1, 0d);
geoInR.AddPoint(dX + 3, dY + 3, 0d);
geoInR.AddPoint(dX + 1, dY + 3, 0d);
geoInR.CloseRings(); //閉じる
//ポリゴン(ポリゴンを構成するリングを追加する)
OGR.Geometry geoPoly = new OGR.Geometry(OGR.wkbGeometryType.wkbPolygon);
geoPoly.AddGeometry(geoOutR);
geoPoly.AddGeometry(geoInR);
//地物を作成→上記のPointを登録&属性値の設定→レイヤに登録
using (OGR.Feature f = new OGR.Feature(fetDef))
{
//属性テーブルのカラムに値を割り当て
f.SetField("fStr", string.Format("ID:{0}", l));
f.SetField("fX", dX);
//地物のジオメトリを設定
f.SetGeometry(geoPoly);
f.SetFID(l);
//レイヤに地物を登録
layer.CreateFeature(f);
}
}
}
}
作図した図はこんな感じで、ちょっと不格好ですが中抜けのポリゴンが出来ているかと思います。
作図した中抜けのポリゴン
2017 Nov. 18追記
ポリゴンを構成するリングを閉じるときは「CloseRings」メソッドを使った方が良さそうです。
(サンプル、55行目、63行目;修正済み)
毎度のオマケ:OGR.Ogr.RegisterAll()で「System.TypeInitializationException」が出たときは、構成マネージャーでx68かx64に変更(64bit版GDALを入れていればx64、32bit版ならx86)すると幸せになれるでしょう。
この記事を作成するにあたって、下記のサイトを参考にさせていただきました。
https://pcjericks.github.io/py-gdalogr-cookbook/geometry.html