7. ジオメトリの操作

ヒント

pyqgisコンソールを使わない場合、このページにあるコードスニペットは次のインポートが必要です:

 1from qgis.core import (
 2  QgsGeometry,
 3  QgsGeometryCollection,
 4  QgsPoint,
 5  QgsPointXY,
 6  QgsWkbTypes,
 7  QgsProject,
 8  QgsFeatureRequest,
 9  QgsVectorLayer,
10  QgsDistanceArea,
11  QgsUnitTypes,
12  QgsCoordinateTransform,
13  QgsCoordinateReferenceSystem
14)

空間地物を表すポイント、ラインストリング、ポリゴンは一般的にジオメトリと呼ばれます。QGISでは、これらは QgsGeometry というクラスで表現されます。

ひとつのジオメトリが実際には単純な(シングルパート)ジオメトリの集合であることがあります。このようなジオメトリは、マルチパートジオメトリと呼ばれます。それが1種類の単純ジオメトリだけを含む場合は、マルチポイント、マルチラインまたはマルチポリゴンと呼びます。例えば、複数の島からなる国は、マルチポリゴンとして表現できます。

ジオメトリの座標値はどの座標参照系(CRS)も利用できます。レイヤーから地物を持ってきたときに、ジオメトリの座標値はレイヤーのCRSのものを持つでしょう。

すべての可能なジオメトリの構成と関係の説明と仕様は、OGC Simple Feature Access Standards で高度な詳細が確認できます。

7.1. ジオメトリの構成

PyQGISでは、ジオメトリを作成するためのいくつかのオプションが用意されています:

  • 座標から

    1gPnt = QgsGeometry.fromPointXY(QgsPointXY(1,1))
    2print(gPnt)
    3gLine = QgsGeometry.fromPolyline([QgsPoint(1, 1), QgsPoint(2, 2)])
    4print(gLine)
    5gPolygon = QgsGeometry.fromPolygonXY([[QgsPointXY(1, 1),
    6    QgsPointXY(2, 2), QgsPointXY(2, 1)]])
    7print(gPolygon)
    

    座標は QgsPoint クラスまたは QgsPointXY クラスを用いて与えられます。これらのクラスの違いは、 QgsPoint はMとZの次元をサポートしていることです。

    ポリライン(Linestring)は、ポイントのリストで表現されます。

    ポリゴンは、線形のリング(すなわち閉じた線分)のリストで表されます。最初のリングは外側のリング(境界)であり、オプションで続くリングはポリゴンの穴となります。いくつかのプログラムとは異なり、QGISはリングを閉じてくれるので、最初のポイントを最後のポイントとして複製する必要はありません。

    マルチパートジオメトリはさらに上のレベルです: マルチポイントはポイントのリストで、マルチラインストリングはラインストリングのリストで、マルチポリゴンはポリゴンのリストです。

  • well-knownテキスト(WKT)から

    geom = QgsGeometry.fromWkt("POINT(3 4)")
    print(geom)
    
  • well-knownバイナリ(WKB)から

    1g = QgsGeometry()
    2wkb = bytes.fromhex("010100000000000000000045400000000000001440")
    3g.fromWkb(wkb)
    4
    5# print WKT representation of the geometry
    6print(g.asWkt())
    

7.2. ジオメトリにアクセス

まず、ジオメトリタイプを調べます。使用するのは wkbType() というメソッドです。これは QgsWkbTypes.Type の列挙から値を返します。

1if gPnt.wkbType() == QgsWkbTypes.Point:
2  print(gPnt.wkbType())
3  # output: 1 for Point
4if gLine.wkbType() == QgsWkbTypes.LineString:
5  print(gLine.wkbType())
6  # output: 2 for LineString
7if gPolygon.wkbType() == QgsWkbTypes.Polygon:
8  print(gPolygon.wkbType())
9  # output: 3 for Polygon

代替案として、 type() メソッドを使うことができます。このメソッドは QgsWkbTypes.GeometryType 列挙から値を返します。

displayString() 関数を使用すると、人間が読めるジオメトリタイプを取得することができます。

1print(QgsWkbTypes.displayString(gPnt.wkbType()))
2# output: 'Point'
3print(QgsWkbTypes.displayString(gLine.wkbType()))
4# output: 'LineString'
5print(QgsWkbTypes.displayString(gPolygon.wkbType()))
6# output: 'Polygon'
Point
LineString
Polygon

また、ジオメトリがマルチパートであるかどうかを調べるヘルパー関数 isMultipart() も用意されています。

ジオメトリから情報を抽出するために、各ベクタ型に対応したアクセサ関数が用意されています。ここでは、これらのアクセサの使い方の例を紹介します:

1print(gPnt.asPoint())
2# output: <QgsPointXY: POINT(1 1)>
3print(gLine.asPolyline())
4# output: [<QgsPointXY: POINT(1 1)>, <QgsPointXY: POINT(2 2)>]
5print(gPolygon.asPolygon())
6# output: [[<QgsPointXY: POINT(1 1)>, <QgsPointXY: POINT(2 2)>, <QgsPointXY: POINT(2 1)>, <QgsPointXY: POINT(1 1)>]]

注釈

タプル (x,y) は実際のタプルではなく、 QgsPoint オブジェクトであり、値は x()y() というメソッドでアクセスできます。

マルチパートジオメトリについては、同様のアクセス関数があります: asMultiPoint(), asMultiPolyline() および asMultiPolygon()

ジオメトリのタイプに関係なく、ジオメトリのすべてのパーツに対して反復処理を行うことが可能です。例

geom = QgsGeometry.fromWkt( 'MultiPoint( 0 0, 1 1, 2 2)' )
for part in geom.parts():
  print(part.asWkt())
Point (0 0)
Point (1 1)
Point (2 2)
geom = QgsGeometry.fromWkt( 'LineString( 0 0, 10 10 )' )
for part in geom.parts():
  print(part.asWkt())
LineString (0 0, 10 10)
gc = QgsGeometryCollection()
gc.fromWkt('GeometryCollection( Point(1 2), Point(11 12), LineString(33 34, 44 45))')
print(gc[1].asWkt())
Point (11 12)

また、 QgsGeometry.parts() メソッドを使ってジオメトリの各パーツを修正することも可能です。

1geom = QgsGeometry.fromWkt( 'MultiPoint( 0 0, 1 1, 2 2)' )
2for part in geom.parts():
3  part.transform(QgsCoordinateTransform(
4    QgsCoordinateReferenceSystem("EPSG:4326"),
5    QgsCoordinateReferenceSystem("EPSG:3111"),
6    QgsProject.instance())
7  )
8
9print(geom.asWkt())
MultiPoint ((-10334728.12541878595948219 -5360106.25905461423099041),(-10462135.16126426123082638 -5217485.4735023295506835),(-10589399.84444035589694977 -5072021.45942386891692877))

7.3. ジオメトリの述語と操作

QGISでは、ジオメトリ述語(contains(), intersects() , ...)や集合演算 (combine(), difference(), ...)などの高度なジオメトリ操作にGEOSライブラリを使っています。 また、面積(ポリゴンの場合)や長さ(ポリゴンとラインの場合)のようなジオメトリのプロパティを計算することもできます。

ここでは、与えられたレイヤーの地物を繰り返し処理し、そのジオメトリに基づいて幾何学的な計算を実行する例を見てみましょう。以下のコードは、チュートリアル QGIS プロジェクトの countries レイヤーにある各国の面積と周囲長を計算し、表示するものです。

以下のコードは レイヤQgsVectorLayer オブジェクトで ポリゴン地物タイプを持っていると仮定します。

 1# let's access the 'countries' layer
 2layer = QgsProject.instance().mapLayersByName('countries')[0]
 3
 4# let's filter for countries that begin with Z, then get their features
 5query = '"name" LIKE \'Z%\''
 6features = layer.getFeatures(QgsFeatureRequest().setFilterExpression(query))
 7
 8# now loop through the features, perform geometry computation and print the results
 9for f in features:
10  geom = f.geometry()
11  name = f.attribute('NAME')
12  print(name)
13  print('Area: ', geom.area())
14  print('Perimeter: ', geom.length())
1Zambia
2Area:  62.822790653431014
3Perimeter:  50.65232014052552
4Zimbabwe
5Area:  33.41113559136511
6Perimeter:  26.608288555013935

さて、これでジオメトリの面積と外周を計算し、印刷することができました。しかし、その値がおかしなことにすぐに気づくでしょう。これは、 QgsGeometry クラスの area()length() を使って計算した場合、面積や周長はCRSを考慮しないためです。より強力な面積と距離の計算を行うには、 QgsDistanceArea クラスを使用することができ、楕円体ベースの計算を行うことができます:

以下のコードは レイヤQgsVectorLayer オブジェクトで ポリゴン地物タイプを持っていると仮定します。

 1d = QgsDistanceArea()
 2d.setEllipsoid('WGS84')
 3
 4layer = QgsProject.instance().mapLayersByName('countries')[0]
 5
 6# let's filter for countries that begin with Z, then get their features
 7query = '"name" LIKE \'Z%\''
 8features = layer.getFeatures(QgsFeatureRequest().setFilterExpression(query))
 9
10for f in features:
11  geom = f.geometry()
12  name = f.attribute('NAME')
13  print(name)
14  print("Perimeter (m):", d.measurePerimeter(geom))
15  print("Area (m2):", d.measureArea(geom))
16
17  # let's calculate and print the area again, but this time in square kilometers
18  print("Area (km2):", d.convertAreaMeasurement(d.measureArea(geom), QgsUnitTypes.AreaSquareKilometers))
1Zambia
2Perimeter (m): 5539361.250294601
3Area (m2): 751989035032.9031
4Area (km2): 751989.0350329031
5Zimbabwe
6Perimeter (m): 2865021.3325076113
7Area (m2): 389267821381.6008
8Area (km2): 389267.8213816008

また、2点間の距離を知りたい場合もあります。

 1d = QgsDistanceArea()
 2d.setEllipsoid('WGS84')
 3
 4# Let's create two points.
 5# Santa claus is a workaholic and needs a summer break,
 6# lets see how far is Tenerife from his home
 7santa = QgsPointXY(25.847899, 66.543456)
 8tenerife = QgsPointXY(-16.5735, 28.0443)
 9
10print("Distance in meters: ", d.measureLine(santa, tenerife))

QGISに含まれているアルゴリズムの多くの例を見つけて、これらのメソッドをベクターデータを分析し変換するために使用できます。ここにそれらのコードのいくつかへのリンクを記載します。