7. ジオメトリの操作

ヒント

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

from qgis.core import (
  QgsGeometry,
  QgsGeometryCollection,
  QgsPoint,
  QgsPointXY,
  QgsWkbTypes,
  QgsProject,
  QgsFeatureRequest,
  QgsVectorLayer,
  QgsDistanceArea,
  QgsUnitTypes,
  QgsCoordinateTransform,
  QgsCoordinateReferenceSystem
)

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

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

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

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

7.1. ジオメトリの構成

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

• 座標から

  gPnt = QgsGeometry.fromPointXY(QgsPointXY(1,1))
  print(gPnt)
  gLine = QgsGeometry.fromPolyline([QgsPoint(1, 1), QgsPoint(2, 2)])
  print(gLine)
  gPolygon = QgsGeometry.fromPolygonXY([[QgsPointXY(1, 1),
      QgsPointXY(2, 2), QgsPointXY(2, 1)]])
  print(gPolygon)
  

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

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

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

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

• well-knownテキスト（WKT）から

  geom = QgsGeometry.fromWkt("POINT(3 4)")
  print(geom)
  

• well-knownバイナリ（WKB）から

  g = QgsGeometry()
  wkb = bytes.fromhex("010100000000000000000045400000000000001440")
  g.fromWkb(wkb)
  
  # print WKT representation of the geometry
  print(g.asWkt())
  

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

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

if gPnt.wkbType() == QgsWkbTypes.Point:
  print(gPnt.wkbType())
  # output: 1 for Point
if gLine.wkbType() == QgsWkbTypes.LineString:
  print(gLine.wkbType())
  # output: 2 for LineString
if gPolygon.wkbType() == QgsWkbTypes.Polygon:
  print(gPolygon.wkbType())
  # output: 3 for Polygon

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

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

print(QgsWkbTypes.displayString(gPnt.wkbType()))
# output: 'Point'
print(QgsWkbTypes.displayString(gLine.wkbType()))
# output: 'LineString'
print(QgsWkbTypes.displayString(gPolygon.wkbType()))
# output: 'Polygon'

Point
LineString
Polygon

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

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

print(gPnt.asPoint())
# output: <QgsPointXY: POINT(1 1)>
print(gLine.asPolyline())
# output: [<QgsPointXY: POINT(1 1)>, <QgsPointXY: POINT(2 2)>]
print(gPolygon.asPolygon())
# 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() メソッドを使ってジオメトリの各パーツを修正することも可能です。

geom = QgsGeometry.fromWkt( 'MultiPoint( 0 0, 1 1, 2 2)' )
for part in geom.parts():
  part.transform(QgsCoordinateTransform(
    QgsCoordinateReferenceSystem("EPSG:4326"),
    QgsCoordinateReferenceSystem("EPSG:3111"),
    QgsProject.instance())
  )

print(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 オブジェクトで ポリゴン地物タイプを持っていると仮定します。

# let's access the 'countries' layer
layer = QgsProject.instance().mapLayersByName('countries')[0]

# let's filter for countries that begin with Z, then get their features
query = '"name" LIKE \'Z%\''
features = layer.getFeatures(QgsFeatureRequest().setFilterExpression(query))

# now loop through the features, perform geometry computation and print the results
for f in features:
  geom = f.geometry()
  name = f.attribute('NAME')
  print(name)
  print('Area: ', geom.area())
  print('Perimeter: ', geom.length())

Zambia
Area:  62.822790653431014
Perimeter:  50.65232014052552
Zimbabwe
Area:  33.41113559136511
Perimeter:  26.608288555013935

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

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

d = QgsDistanceArea()
d.setEllipsoid('WGS84')

layer = QgsProject.instance().mapLayersByName('countries')[0]

# let's filter for countries that begin with Z, then get their features
query = '"name" LIKE \'Z%\''
features = layer.getFeatures(QgsFeatureRequest().setFilterExpression(query))

for f in features:
  geom = f.geometry()
  name = f.attribute('NAME')
  print(name)
  print("Perimeter (m):", d.measurePerimeter(geom))
  print("Area (m2):", d.measureArea(geom))

  # let's calculate and print the area again, but this time in square kilometers
  print("Area (km2):", d.convertAreaMeasurement(d.measureArea(geom), QgsUnitTypes.AreaSquareKilometers))

Zambia
Perimeter (m): 5539361.250294601
Area (m2): 751989035032.9031
Area (km2): 751989.0350329031
Zimbabwe
Perimeter (m): 2865021.3325076113
Area (m2): 389267821381.6008
Area (km2): 389267.8213816008

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

d = QgsDistanceArea()
d.setEllipsoid('WGS84')

# Let's create two points.
# Santa claus is a workaholic and needs a summer break,
# lets see how far is Tenerife from his home
santa = QgsPointXY(25.847899, 66.543456)
tenerife = QgsPointXY(-16.5735, 28.0443)

print("Distance in meters: ", d.measureLine(santa, tenerife))

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

• QgsDistanceArea クラス: Distance matrix algorithm による距離と領域の計算

• Lines to polygons algorithm