PyQGIS ラスタから値を取得するサンプル


どんな動きをするの?

f:id:Chiakikun:20200617221955p:plain

ラスタレイヤをアクティブにした状態でプラグインを実行すると、マウスカーソル位置の標高をPythonコンソールに表示します。


コード

import qgis
from qgis.core import *
from qgis.gui  import *
import pyproj

class GetRasterPixelValue:

    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)
    

ソースはこちら


使い方

例として、マウスイベントサンプルに追加します。このmouseeventsample.pyの使い方はこちらを参照してください。

chiakikun.hatenadiary.com

mouseeventsample.pyとgetrasterpixelvalue.pyを同じフォルダに置いて、mouseeventsample.pyのインポート部分に次のコードを追記します。

from .getrasterpixelvalue import GetRasterPixelValue

mouseeventsample.pyの__init__の引数は次のように変更してください。

__init__(self, iface, canvas)

mouseeventsample.pyの__init__に次のコードを加えます。

self.gr = GetRasterPixelValue(iface.activeLayer())

mouseeventsample.pyのcanvasMoveEventを次のコードに変えます。

def canvasMoveEvent(self, event):
    print(str(self.gr.getValueInterpolation(self.toMapCoordinates(event.pos()))))

nodialog_skelton.pyのstartを次のコードに書き換えます。

    def start(self):
        maptool = MouseEventSample(self.iface, self.canvas)

nodialog_skelton.pyはこのようになればOKです。

~略~

from .mouseeventsample import MouseEventSample

class NodialogSkelton(qgis.gui.QgsMapTool):

    def start(self):
        maptool = MouseEventSample(self.iface, self.canvas)

        self.canvas.setMapTool(maptool)
        self.canvas.mapToolSet.connect(self.unsetTool) # このサンプル実行中に他のアイコンを押した場合
~略~

 

mouseeventsample.pyはこのようになればOKです。

~略~
from .getrasterpixelvalue import GetRasterPixelValue

class MouseEventSample(QgsMapTool):

    def __init__(self, iface, canvas):
        QgsMapTool.__init__(self, canvas)

        self.gr = GetRasterPixelValue(iface.activeLayer())


    def canvasMoveEvent(self, event):
        print(str(self.gr.getValueInterpolation(self.toMapCoordinates(event.pos()))))
~略~

最後までご覧頂き、ありがとうございました