QGIS_24_directions.py

# -*- coding: utf-8 -*-


import math
from qgis.PyQt.QtCore import (QCoreApplication, QVariant)
from qgis.core import (QgsPointXY, QgsGeometry, QgsFeature,
        QgsField, QgsFields, QgsWkbTypes, QgsProject, QgsUnitTypes,
        QgsCoordinateReferenceSystem, QgsCoordinateTransform)
from qgis.core import (QgsProcessing,
                       QgsFeatureSink,
                       QgsProcessingException,
                       QgsProcessingAlgorithm,
                       QgsProcessingParameterDefinition,
                       QgsProcessingParameterPoint,
                       QgsProcessingParameterDistance,
                       QgsProcessingParameterFeatureSink)
from qgis import processing


class Qgis24DirectionsAlgorithm(QgsProcessingAlgorithm):
    """
    """
    
    # 各パラメータの ID
    POINT = 'POINT'
    DISTANCE = 'DISTANCE'
    OUTPUT = 'OUTPUT'
    
    DIRECTIONS = r'子癸丑艮寅甲卯乙辰巽巳丙午丁未坤申庚酉辛戌乾亥壬'
    
    def tr(self, string):
        return QCoreApplication.translate('Processing', string)
    
    def createInstance(self):
        return Qgis24DirectionsAlgorithm()
    
    def name(self):
        # アルゴリズムの ID
        return '24 directions'
    
    def displayName(self):
        # アルゴリズムが表示されるときの名称
        return self.tr('24 directions')
    
    def shortHelpString(self):
        # アルゴリズムのヘルプとして表示される文章
        return self.tr("generate 24 directions")
    
    def initAlgorithm(self, config=None):
        # パラメータの定義を行う
        
        self.addParameter(
            QgsProcessingParameterPoint(
                self.POINT,
                self.tr('Start point')
            )
        )
        
        dist_param = QgsProcessingParameterDistance(
            self.DISTANCE,
            self.tr('Distance'),
            defaultValue=100, # km
            # パラメータ定義時の minValue は単位に関わらず適用されるため
            # 最低値のチェック及び適用はあとで行う
            # とはいえウィジットレベルで最低0以上の制約は行う
            minValue=0
        )
        dist_param.setDefaultUnit(QgsUnitTypes.DistanceKilometers)
        # 高度なパラメータフラグ
        dist_param.setFlags(QgsProcessingParameterDefinition.FlagAdvanced)
        self.addParameter(dist_param)
        
        # 結果の出力
        self.addParameter(
            QgsProcessingParameterFeatureSink(
                self.OUTPUT,
                self.tr('Output layer')
            )
        )
    
    def processAlgorithm(self, parameters, context, feedback):
        # 実行される処理を記述する
        
        # self.parameterAsFoo でパラメータを参照できる
        # 詳しくは QgsProcessingAlgorithm クラスのドキュメント参照
        
        crs4326 = QgsCoordinateReferenceSystem('EPSG:4326')
        
        # パラメータの値を参照する（座標値は EPSG:4326 に変換される）
        point = self.parameterAsPoint(
            parameters,
            self.POINT,
            context,
            crs4326
        )
        if point is None:
            raise QgsProcessingException(
                'Error: failed to get coordnates of POINT')
        # 実行時のログ欄に出力
        feedback.pushInfo('point ({}, {})'.format(point.x(), point.y()))
        
        # QgsProcessingParameterDistance の場合は、デフォルトの単位に
        # 変換された数値が得られる？
        dist = self.parameterAsDouble(
            parameters,
            self.DISTANCE,
            context
        )
        # メートルになおす
        # addParameter 時の minValue は単位に関わらず適用されるので
        # 最低値等をここで設定する
        dist = max((10, dist)) * 1000.0
        feedback.pushInfo('dist: {}'.format(dist))
        
        # 出力レイヤの属性の定義を行う
        fields = QgsFields()
        fields.append(QgsField('id', QVariant.Int))
        fields.append(QgsField('lat', QVariant.Double))
        fields.append(QgsField('lng', QVariant.Double))
        fields.append(QgsField('direction', QVariant.String))
        
        (sink, dest_id) = self.parameterAsSink(
            parameters,
            self.OUTPUT,
            context,
            fields,
            QgsWkbTypes.Polygon,
            crs4326
        )
        
        if sink is None:
            raise QgsProcessingException(self.invalidSinkError(parameters, self.OUTPUT))
        
        
        # 点を中心とした正距方位図法の投影法を定義する
        crs_aeqd = QgsCoordinateReferenceSystem()
        crs_aeqd.createFromProj(
            '+proj=aeqd +lat_0={} +lon_0={} +datum=WGS84'.format(
                point.y(), point.x()))
        if not crs_aeqd.isValid():
            raise QgsProcessingException('Invalid CRS')
        
        trans = QgsCoordinateTransform(crs4326, crs_aeqd,
            QgsProject.instance())
        if not trans.isValid():
            raise QgsProcessingException('Invalid coordinate transformer')
        
        
        # セグメントに分割することで曲線のような折れ線を作る
        segment_dist = 10000
        num_segments = dist / segment_dist
        if num_segments - math.floor(num_segments) > 0:
            num_segments += 1
        num_segments = int(num_segments)
        
        for theta in range(0, 24):
            t = (6 - theta) % 24
            # ポリゴン上の点
            points = []
            
            cost = math.cos((t - 0.5) * 15 / 180 * math.pi)
            sint = math.sin((t - 0.5) * 15 / 180 * math.pi)
            for i in range(0, num_segments):
                # 点中心正距方位図法における、線上の点を
                # EPSG:4326 に逆変換を行う
                x = i * segment_dist * cost
                y = i * segment_dist * sint
                p_4326 = trans.transform(QgsPointXY(x, y),
                    QgsCoordinateTransform.ReverseTransform)
                points.append(p_4326)
            p_4326 = trans.transform(QgsPointXY(dist*cost, dist*sint),
                QgsCoordinateTransform.ReverseTransform)
            points.append(p_4326)
            
            
            cost = math.cos((t + 0.5) * 15 / 180 * math.pi)
            sint = math.sin((t + 0.5) * 15 / 180 * math.pi)
            p_4326 = trans.transform(QgsPointXY(dist*cost, dist*sint),
                QgsCoordinateTransform.ReverseTransform)
            points.append(p_4326)
            for i in range(0, num_segments):
                # 点中心正距方位図法における、線上の点を
                # EPSG:4326 に逆変換を行う
                x = (num_segments - i) * segment_dist * cost
                y = (num_segments - i) * segment_dist * sint
                p_4326 = trans.transform(QgsPointXY(x, y),
                    QgsCoordinateTransform.ReverseTransform)
                points.append(p_4326)
            points.append(points[0])
            
            
            # 折れ線の測地線のジオメトリを作る
            geom = QgsGeometry().fromPolygonXY([points])
            
            # 折れ線の測地線のジオメトリを持つ地物を作る
            feat = QgsFeature(fields)
            feat.setGeometry(geom)
            feat.setAttribute('id', theta)
            feat.setAttribute('lat', point.y())
            feat.setAttribute('lng', point.x())
            feat.setAttribute('direction', self.DIRECTIONS[theta])
            
            # 地物を出力レイヤ (sink) に追加する
            sink.addFeature(feat, QgsFeatureSink.FastInsert)
        
        # アルゴリズムの結果を返す
        # このアルゴリズムでは出力 sink のひとつのみ
        return {self.OUTPUT: dest_id}