PostGIS ST_Intersectsで地物を抽出する
はじめに
PostGISを使って、北海道に重なる基準地域メッシュ（第３次地域区画）のオブジェクト抽出を行ったので、その時の作業メモです。
メッシュ（PostGISに「jpmesh_3」で登録、フィールド「code」にはメッシュ番号が入っています）
拡大したところ。
北海道（PostGISに「北海道」で登録）
コード
INSERT INTO
temp1(geom, メッシュコード)
SELECT
jpmesh_3.geom,
jpmesh_3.code
FROM
jpmesh_3, 北海道
WHERE
ST_Intersects(
jpmesh_3.geom,北海道.geom
);
上記のコードで、北海道に重なるメッシュのオブジェクトを抽出して、予め作成しておいたtemp1というテーブルに投入しています。
実行速度を計測
ネタのかさ増しのため、INDEX作成前と後で実行速度がどのくらい変わるのか計測してみました。
INDEX作成前 Time: 6547.685 ms (00:06.548)
INDEX作成後 Time: 4210.643 ms (00:04.211)
ついでに北海道のオブジェクトのノードを半分程間引いて計測してみました。
間引いた後の北海道はこんな感じです
↓
↓
赤く見えているところが、間引き前の北海道です。
INDEX作成前 Time: 4933.373 ms (00:04.933)
INDEX作成後 Time: 3123.348 ms (00:03.123)
QGISで作業する場合
QGISで同じことをする場合は「場所による選択」を実行します。
実行時間は121.43秒でした。PostGISにテーブルを登録する手間を考えると、QGISでの作業は実行時間こそかかりましたが、簡単操作でできるので、少ないレコード数のレイヤを抽出するならQGISで...ということになってしまいました。
最後までご覧いただき、ありがとうございました。
QGISのプラグイン作成 オブジェクトを移動、拡大縮小、回転させるサンプル
ここで紹介しているサンプルはオブジェクトの形状を変更するものなので、実行する際はバックアップを必ず取ってください。
はじめに
通常、マルチオブジェクトを移動したり回転させようとするとパーツ全体が影響をうけますが、このサンプルでは1つのパーツだけを移動、拡大縮小、回転させることができます。
使い方
ここからzipファイルをダウンロードしてインストールします。インストールできたら、処理内容に合わせてmultiobjectedit_sample.pyの以下の部分を編集してください。
self.maptool = MultiObjectEditClass(self.iface, self.canvas, 'move') # 移動の場合
3番目の引数が、移動させたいなら'move'、拡大縮小させたいなら'scale'、回転させたいなら'rotate'に設定してください。
# self.layer.startEditing() # バックアップは取りましたか？バックアップできたなら、#を外してください
self.changeGeometry(self.layer.selectedFeatures()[0])
# self.layer.commitChanges() # ここも
レイヤをバックアップしたら、#を外してください。このままだと編集内容が反映されません。
レイヤをアクティブにして、プラグインを実行します。オブジェクトをクリックすると...
オブジェクトが編集状態になり、移動、拡大縮小、回転ができるようになります。上の画像では、マルチポリゴンの一部を回転させています。
回転できました。オブジェクトから少し離れたところで右クリックすると、編集状態を解除できます。
ライン（シングルオブジェクト）を拡大しています。
マルチポリゴンの一部を移動させています。
注意 このサンプルではドーナツのポリゴンには対応していません。
コード
class MultiObjectEditSample(QgsMapTool):
def selectedFeature(self, feature):
self.maptool.setFeature(feature)
def start(self):
self.maptool = MultiObjectEditClass(self.iface, self.canvas, 'move') # 移動の場合
self.maptool.setLayer(self.iface.activeLayer())
self.maptool.featureIdentified.connect(self.selectedFeature)
～略～
class MultiObjectEditClass(QgsMapToolIdentifyFeature):
setFeature = pyqtSignal(QgsGeometry)
def __init__(self, iface, canvas, edittype): # edittype = 'scale', 'move', 'rotate'
self.iface = iface
self.canvas = canvas
self.layer = self.iface.activeLayer()
self.objType = self.layer.geometryType()
self.edittype = edittype
QgsMapToolIdentifyFeature.__init__(self, self.canvas)
self.myRubberBand = None
def canvasPressEvent(self, event): # featureIdentifiedより先に呼ばれる
self.srcpos = event.pos()
if self.myRubberBand == None: return
if event.button() == QtCore.Qt.LeftButton:
# self.layer.startEditing() # バックアップは取りましたか？バックアップできたなら、#を外してください
self.changeGeometry(self.layer.selectedFeatures()[0])
# self.layer.commitChanges() # ここも
self.myRubberBand.reset()
self.myRubberBand = None
self.layer.removeSelection()
def setFeature(self, feat):
if self.myRubberBand != None: return
self.layer.select(feat.id())
curpos = QgsGeometry().fromPointXY(self.toMapCoordinates(self.srcpos))
mpol = self.tolist(feat.geometry())
# 選択したマルチオブジェクトで、選択位置に一番近いオブジェクト取得
dist = [pol.distance(curpos) for pol in mpol]
self.nearidx = dist.index(min(dist))
self.nearobj = mpol[self.nearidx]
# ラバーバンドにnearobjを追加する
self.myRubberBand = QgsRubberBand( self.canvas, self.objType )
self.myRubberBand.setColor(QColor(255, 0, 0, 255))
self.myRubberBand.addGeometry(self.nearobj, self.layer)
def canvasMoveEvent(self, e):
if self.myRubberBand == None: return
self.myRubberBand.reset(self.objType)
dstpos = e.pos()
if self.edittype == 'move':
self.moveObject(self.srcpos, dstpos)
elif self.edittype == 'scale':
self.scaleObject(self.srcpos, dstpos)
else:
self.rotateObject(self.srcpos, dstpos)
def changeGeometry(self, feat):
if self.objType == QgsWkbTypes.GeometryType.PolygonGeometry:
if feat.geometry().isMultipart():
mpol = feat.geometry().asMultiPolygon()
mpol[self.nearidx] = self.myRubberBand.asGeometry().asPolygon()
self.layer.changeGeometry (feat.id(), QgsGeometry().fromMultiPolygonXY(mpol))
else:
pol = self.myRubberBand.asGeometry().asPolygon()
self.layer.changeGeometry (feat.id(), QgsGeometry.fromPolygonXY(pol))
elif self.objType == QgsWkbTypes.GeometryType.LineGeometry:
if feat.geometry().isMultipart():
mline = feat.geometry().asMultiPolyline()
mline[self.nearidx] = self.myRubberBand.asGeometry().asPolyline()
self.layer.changeGeometry (feat.id(), QgsGeometry().fromMultiPolylineXY(mline))
else:
line = self.myRubberBand.asGeometry().asPolyline()
self.layer.changeGeometry (feat.id(), QgsGeometry.fromPolylineXY(line))
elif self.objType == QgsWkbTypes.GeometryType.PointGeometry:
if feat.geometry().isMultipart():
mpnt = feat.geometry().asMultiPoint()
mpnt[self.nearidx] = self.myRubberBand.asGeometry().asMultiPoint()[0]
self.layer.changeGeometry (feat.id(), QgsGeometry().fromMultiPointXY(mpnt))
else:
pnt = self.myRubberBand.asGeometry().asMultiPoint()[0]
self.layer.changeGeometry (feat.id(), QgsGeometry.fromPointXY(pnt))
def tolist(self, obj):
if self.objType == QgsWkbTypes.GeometryType.PolygonGeometry:
if obj.isMultipart():
return [QgsGeometry().fromPolygonXY(pol) for pol in obj.asMultiPolygon()]
else:
return [QgsGeometry().fromPolygonXY(obj.asPolygon())]
elif self.objType == QgsWkbTypes.GeometryType.LineGeometry:
if obj.isMultipart():
return [QgsGeometry().fromPolylineXY(pline) for pline in obj.asMultiPolyline()]
else:
return [QgsGeometry().fromPolylineXY(obj.asPolyline())]
elif self.objType == QgsWkbTypes.GeometryType.PointGeometry:
if obj.isMultipart():
return [QgsGeometry().fromPointXY(pnt) for pnt in obj.asMultiPoint()]
else:
return [QgsGeometry().fromPointXY(obj.asPoint())]
def moveObject(self, srcpos, dstpos):
temp = QgsGeometry(self.nearobj)
s = self.toMapCoordinates(srcpos)
d = self.toMapCoordinates(dstpos)
self.myRubberBand.addGeometry(self.trans(temp, s, d), self.layer)
def scaleObject(self, srcpos, dstpos):
original = QgsPointXY(0, 0)
centroid = self.nearobj.centroid().asPoint()
temp = QgsGeometry(self.nearobj)
s = QgsGeometry().fromPointXY(QgsPointXY(srcpos))
d = QgsGeometry().fromPointXY(QgsPointXY(dstpos))
# 一度中心に持ってこないとオブジェクトが変な方向に行ってしまうので
temp = self.trans(temp, centroid, original)
temp = self.scale(temp, s, d)
temp = self.trans(temp, original, centroid)
self.myRubberBand.addGeometry(temp, self.layer)
def rotateObject(self, srcpos, dstpos):
original = QgsPointXY(0, 0)
centroid = self.nearobj.centroid().asPoint()
temp = QgsGeometry(self.nearobj)
s = QgsPointXY(self.toMapCoordinates(srcpos))
d = QgsPointXY(self.toMapCoordinates(dstpos))
# 一度中心に持ってこないとオブジェクトが変な方向に行ってしまうので
temp = self.trans(temp, centroid, original)
temp = self.rotate(temp, s, d)
temp = self.trans(temp, original, centroid)
self.myRubberBand.addGeometry(temp, self.layer)
def trans(self, obj, srcpos, dstpos):
x = dstpos.x() - srcpos.x()
y = dstpos.y() - srcpos.y()
obj.transform(QTransform(1, 0, 0, 1, x, y))
return obj
def scale(self, obj, srcpos, dstpos):
cent = QgsGeometry().fromPointXY(QgsPointXY(0, 0))
srclen = cent.distance(srcpos)
dstlen = cent.distance(dstpos)
ratio = (dstlen / srclen) * 1.5 # 「1.5」は適当
obj.transform(QTransform(ratio, 0, 0, ratio, 0, 0))
return obj
def rotate(self, obj, srcpos, dstpos):
centroid = self.nearobj.centroid().asPoint()
theta = self.radian(centroid, dstpos)
sin = math.sin(theta)
cos = math.cos(theta)
obj.transform(QTransform(cos, sin, -sin, cos, 0, 0))
return obj
def radian(self, a, b):
return math.atan2(b.y() - a.y(), b.x() - a.x())
def deactivate(self):
try: # 左クリックで既にNoneしていた場合エラーになるので
self.myRubberBand.reset()
except:
pass
self.myRubberBand = None
self.layer.removeSelection()
8行目 編集処理内容を設定します。
9行目 編集するレイヤを設定します。
4行目 地物を選択した時に、編集を開始します。
50～55行目 選択したマルチオブジェクトをバラして、マウスカーソルに一番近いパーツを編集対象とします。
70～75行目 マウスカーソルを動かして、移動、拡大縮小、回転をさせています。
33行目 左ボタンで編集内容を確定します。
コード全文
最後までご覧いただき、ありがとうございました。
QGISのプラグイン作成 ラスタから値を補間して取得するサンプル
はじめに
このプラグインのサンプルでは、DEMのラスタから標高値を補完して取得するサンプルです。
使い方
ここからzipファイルをダウンロードしてインストールします。
ラスタレイヤをアクティブにした状態でプラグインを実行すると、マウスカーソル位置の標高をPythonコンソールに表示します。
コード
class GetRasterPixelValueSample(QgsMapTool):
def canvasMoveEvent(self, event):
print(str(self.gr.getValueInterpolation(self.toMapCoordinates(event.pos()))))
def start(self):
self.gr = GetRasterPixelValueClass(self.iface.activeLayer())
～略～
class GetRasterPixelValueClass:
def __init__(self, layer: QgsRasterLayer):
self.__grs80 = pyproj.Geod(ellps='GRS80')
# ピクセル1つのサイズ
self.__unitX = layer.rasterUnitsPerPixelX()
self.__unitY = layer.rasterUnitsPerPixelY()
# ラスタのサイズ
self.__height = layer.height()
self.__width = layer.width()
# ラスタの原点
self.__provider = layer.dataProvider()
self.__extent = self.__provider.extent()
self.__orgX = self.__extent.xMinimum()
self.__orgY = self.__extent.yMaximum()
# ピクセル
self.block = self.__provider.block(1, self.__extent, self.__width, self.__height)
def getIndex(self, pos: QgsPointXY):
pixelIdxX = (pos.x() - self.__orgX) / self.__unitX
pixelIdxY = (self.__orgY - pos.y()) / self.__unitY
if (pixelIdxX < 0) or (pixelIdxY < 0):
return None
pixelIdxX = int(pixelIdxX)
pixelIdxY = int(pixelIdxY)
if (pixelIdxX > self.__width - 1) or (pixelIdxY > self.__height - 1):
return None
return [pixelIdxX, pixelIdxY]
def getGeometry(self, idxX: int, idxY: int):
if (idxX < 0) or (idxY < 0):
return None
if (idxX > self.__width - 1) or (idxY > self.__height - 1):
return None
pixelPosX = self.__unitX * idxX + self.__orgX + self.__unitX / 2
pixelPosY = self.__orgY - self.__unitY * idxY - self.__unitY / 2
return [pixelPosX, pixelPosY]
def getValueFromPos(self, pos: QgsPointXY):
ident = self.__provider.identify(pos, QgsRaster.IdentifyFormatValue )
if ident.isValid():
return list(ident.results().values())[0]
else:
return None
def getValueFromIdx(self, idxX: int, idxY: int):
if (idxX < 0) or (idxY < 0):
return None
if (idxX > self.__width - 1) or (idxY > self.__height - 1):
return None
return self.block.value(idxY, idxX)
def getValueInterpolation(self, pos: QgsPointXY):
pixidx = self.getIndex(pos)
if pixidx == None:
return None
pixelIdxX = pixidx[0]
pixelIdxY = pixidx[1]
pixpos = self.getGeometry(pixelIdxX, pixelIdxY)
pixelPosX = pixpos[0]
pixelPosY = pixpos[1]
# 右上
if (pixelPosX < pos.x()) and (pixelPosY < pos.y()):
if (pixelIdxX > self.__width - 2) or (pixelIdxY - 1 < 0):
return None
a = self.__getXYZ__(pixelIdxX+1, pixelIdxY-1)
b = self.__getXYZ__(pixelIdxX, pixelIdxY-1)
c = self.__getXYZ__(pixelIdxX, pixelIdxY)
d = self.__getXYZ__(pixelIdxX+1, pixelIdxY)
# 左上
elif (pixelPosX >= pos.x()) and (pixelPosY < pos.y()):
if (pixelIdxX - 1 < 0) or (pixelIdxY - 1 < 0):
return None
a = self.__getXYZ__(pixelIdxX, pixelIdxY-1)
b = self.__getXYZ__(pixelIdxX-1, pixelIdxY-1)
c = self.__getXYZ__(pixelIdxX-1, pixelIdxY)
d = self.__getXYZ__(pixelIdxX, pixelIdxY)
# 左下
elif (pixelPosX >= pos.x()) and (pixelPosY >= pos.y()):
if (pixelIdxX - 1 < 0) or (pixelIdxY > self.__height - 2):
return None
a = self.__getXYZ__(pixelIdxX, pixelIdxY)
b = self.__getXYZ__(pixelIdxX-1, pixelIdxY)
c = self.__getXYZ__(pixelIdxX-1, pixelIdxY+1)
d = self.__getXYZ__(pixelIdxX, pixelIdxY+1)
# 右下
elif (pixelPosX < pos.x()) and (pixelPosY >= pos.y()):
if (pixelIdxX > self.__width - 2) or (pixelIdxY > self.__height - 2):
return None
a = self.__getXYZ__(pixelIdxX+1, pixelIdxY)
b = self.__getXYZ__(pixelIdxX, pixelIdxY)
c = self.__getXYZ__(pixelIdxX, pixelIdxY+1)
d = self.__getXYZ__(pixelIdxX+1, pixelIdxY+1)
dist1 = self.__grs80.inv(b.x(), b.y(), a.x(), a.y())[2]
dist2 = self.__grs80.inv(b.x(), b.y(), pos.x(), b.y())[2]
delta = (a.z() - b.z()) / dist1
tmp1 = delta * dist2 + b.z()
dist1 = self.__grs80.inv(c.x(), c.y(), d.x(), d.y())[2]
dist2 = self.__grs80.inv(c.x(), c.y(), pos.x(), c.y())[2]
delta = (d.z() - c.z()) / dist1
tmp2 = delta * dist2 + c.z()
dist1 = self.__grs80.inv(d.x(), d.y(), a.x(), a.y())[2]
dist2 = self.__grs80.inv(d.x(), d.y(), d.x(), pos.y())[2]
delta = (tmp1 - tmp2) / dist1
inp = delta * dist2 + tmp2
return inp
def __getXYZ__(self, idxX: int, idxY: int):
x = self.__unitX * idxX + self.__orgX + self.__unitX / 2
y = self.__orgY - self.__unitY * idxY - self.__unitY / 2
z = self.block.value(idxY, idxX)
return QgsPoint(x, y, z)
90～121行目 現在の位置の周囲4ピクセルの座標を決めています。コメントで「右上」とある箇所では、現在のピクセルの上、右上、右、現在位置の4ピクセルを取得します。
145行目 その位置のピクセル値をそれぞれ取得します。
123～136行目 補間して、マウス位置の標高値を取得します。
コード全文
最後までご覧いただき、ありがとうございました。
