ポリゴン内の高さ抽出の件

[GTT]

お世話になっております。

【指定条件の高さ抽出】の件で、ご教授賜りたく連絡させていただきました。

ポリゴンと重畳する点群データから、以下2つの高さをポリゴンの属性に付与したいと考えております。

• ポリゴン内にある点群データの高さの中央値
• ポリゴン内にある点群データの95％の高さ
  • 点が200点ある場合は、高さを昇順に並べ190番目の高さを抽出する、というものでございます。

上記の処理を行うとなると、pythonでの処理になると存じますが、私自身、詳しくございません。
お手数をおかけいたしますが、ご教示賜りますようお願い申し上げます。

最終的には以下の図のように、複数ポリゴンに重畳する点群データから、

【中央値】と【95％高さ】を抽出、ポリゴンに属性付与できればと考えております。

もし試行して賜れる方がいらっしゃるならば、本メールにサンプルデータを添付させていただきますので、ご活用いただけましたら幸いでございます。

[ありた]

こんにちは

任意値のパーセンタイルは対応していませんが、中央値や四分位数であれば QGIS の機能で可能です。

中央値や四分位数を QGIS を使って求めるには、
LAS データは近年 QGIS でも読み込めますが、伝統的な GIS のファイル形式とは異なる
点群形式のデータであるため、解析処理を行うにあたり形式変換が必要です。
また、座標値として XYZ 座標は有していますが、集計のための Z 属性を持っていません。
このため、まず
(1) LAS データと建物ポリゴンデータを QGIS に読み込ませる
(2) LAS データをエクスポートで点ベクタ形式に変換する
(3) プロセシングツールの「Z 値を抽出」で、座標値の Z 座標を属性として付与

データの変換等の整備ができましたら、ポリゴン領域で集計し統計情報の取得処理をします。
(4) 「空間結合（集計つき）」で地物を結合するレイヤに建物ポリゴン、比較対象レイヤを
　Z 属性付与した点群レイヤを指定。集計する属性に Z 属性を指定し、集計関数に中央値を設定
以上で可能かと思います。



パーセンタイルに対応したい場合は、 Python あるいは R などで処理を記述する必要があります。
QGIS 上でも Python を動かすことができますが、読み込んだレイヤに対する処理になるため
個人的には QGIS 上ではなく、素の Python スクリプトで処理した方がわかりやすいと思います。
（ただし R が利用可能なら Processing R Provider を利用する方法もあるかと思います。）

Python および、ある程度のライブラリは QGIS と一緒にインストールされているので、
Windows の場合、スタートメニューのすべてのプログラムより QGIS 3.xx.y フォルダ内にある

OSGeo4W Shell を起動し呼び出すことで可能です。

また前準備として LAS データを Python で読み込む必要があるため、 OSGeo4W Shell 上で

下記コマンドにより laspy ライブラリをインストールします。
python -m pip install -U pip
pip install -U laspy

そのうえ、たとえば下記のような処理をしてやることで可能かと思います。

from shapely.geometry import Point
import numpy as np
import pandas as pd
import geopandas as gpd
import laspy

# 座標参照系の設定は、システムとしては LAS データの設定値を参照する方が妥当だが、
# 特定地域を想定してハードコーディングで対応
# また3次元データとして、高度も含めた EPSG:10169 あるいは EPSG:6676+6695 の方が
# 本質的には適切かもしれないが互換性と、あとでポリゴンと空間解析を行うため
# EPSG:6676 が無難か
#crs = 'EPSG:10169'
#crs = 'EPSG:6676+6695'
crs = 'EPSG:6676'

las = laspy.read("./sample_pointcloud.las")

# LAS データ内の点群のデータから GeoDataFrame を作成
lasdf = gpd.GeoDataFrame(las.points.array)
# X, Y, Z の各値は scale と offset で補正しないと実座標にならない
# https://qiita.com/nokonoko_1203/items/8b8568545148acfe11d3
lasdf['X'] = lasdf['X'] * las.header.x_scale + las.header.x_offset
lasdf['Y'] = lasdf['Y'] * las.header.y_scale + las.header.y_offset
lasdf['Z'] = lasdf['Z'] * las.header.z_scale + las.header.z_offset

# この時点では X, Y, Z は属性データなのでジオメトリデータを作成
lasdf['geom'] = lasdf.apply(lambda r: Point(r['X'], r['Y'], r['Z']), axis=1)
lasdf.set_geometry('geom', inplace=True, crs=crs)


# ポリゴンデータを読み込む
polygon = gpd.read_file("./sample_bldg.shp")
# 念のため、再投影を行い座標参照系を揃える
polygon.to_crs(crs=crs, inplace=True)

# もし既に各ポリゴンを識別できる属性があれば省略してよいが、ない場合は
# 下記により行番号を付与する
polygon['LINENO'] = range(1, len(polygon)+1)


# 建物ポリゴンと点群を空間結合する
# この時点で重なっている点群の数だけ結合されたデータがあるので
# groupby で建物ポリゴンごとにグルーピングして、そのグループで
# Z 属性 (Z 座標) について集計を行う。ここでは中央値と 95% パーセンタイル
agg_result = gpd.sjoin(polygon, lasdf).groupby('LINENO', as_index=False)['Z'].agg(
  {'median': 'median', '95%': lambda x: np.quantile(x, 0.95)})
# 得られた統計量を元の建物ポリゴンに結合させる
polygon = pd.merge(polygon, agg_result, on='LINENO')

# ファイルに保存
polygon.to_file("./bldg_with_stat.gpkg", driver="GPKG")

なお、下記コマンドであらかじめ LAS データを点データに変換を行っている場合には

laspy のインストールおよび上記スクリプトの前半部分（ LAS データに対する処理）は不要になります。

（近年の QGIS であれば pdal は同梱されインストールされているはずです。）

pdal translate -w writers.ogr --writers.ogr.ogrdriver=GPKG --writers.ogr.attr_dims=Z sample_pointcloud.las sample_pointcloud.gpkg

2024年5月20日月曜日 11:30:32 UTC+9 GTT:

[GTT]

ありた様

いつも丁寧なご対応をいただき、誠にありがとうございます。
この度は長文にわたるご説明をさせていただき、本当に助かります。

記載させていただいた内容に沿って対応させていただきます。

取り急ぎ、感謝申し上げます。

2024年5月22日水曜日 1:01:15 UTC+9 ありた:

[GTT]

ありた様

お世話になっております。

ご教示いただいた方法によりまして、無事ほしい結果を抽出することができました。

いつも誠にありがとうございます。

これからも、どうぞよろしくお願い申し上げます。

2024年5月22日水曜日 11:19:13 UTC+9 GTT:

[GTT]

ありた様

いつもお世話になっております。

【OSgeo4W Shell上で行うlaspyファイルのインストール】について、連絡させて頂きました。

いちどご教授頂いた処理方法にて、希望していた出力結果を出すことが出来たのですが、

日を開けて再度処理をしようと、OSgeo4W Shell上でlaspyファイルのインストール行うとすると以下のエラーが起きました。

C:\Program Files\QGIS 3.34.5>pip install -U laspy
Defaulting to user installation because normal site-packages is not writeable
WARNING: Retrying (Retry(total=4, connect=None, read=None, redirect=None, status=None)) after connection broken by 'SSLError(SSLCertVerificationError(1, '[SSL: CERTIFICATE_VERIFY_FAILED] certificate verify failed: unable to get local issuer certificate (_ssl.c:1129)'))': /simple/laspy/
WARNING: Retrying (Retry(total=3, connect=None, read=None, redirect=None, status=None)) after connection broken by 'SSLError(SSLCertVerificationError(1, '[SSL: CERTIFICATE_VERIFY_FAILED] certificate verify failed: unable to get local issuer certificate (_ssl.c:1129)'))': /simple/laspy/
WARNING: Retrying (Retry(total=2, connect=None, read=None, redirect=None, status=None)) after connection broken by 'SSLError(SSLCertVerificationError(1, '[SSL: CERTIFICATE_VERIFY_FAILED] certificate verify failed: unable to get local issuer certificate (_ssl.c:1129)'))': /simple/laspy/
WARNING: Retrying (Retry(total=1, connect=None, read=None, redirect=None, status=None)) after connection broken by 'SSLError(SSLCertVerificationError(1, '[SSL: CERTIFICATE_VERIFY_FAILED] certificate verify failed: unable to get local issuer certificate (_ssl.c:1129)'))': /simple/laspy/
WARNING: Retrying (Retry(total=0, connect=None, read=None, redirect=None, status=None)) after connection broken by 'SSLError(SSLCertVerificationError(1, '[SSL: CERTIFICATE_VERIFY_FAILED] certificate verify failed: unable to get local issuer certificate (_ssl.c:1129)'))': /simple/laspy/
Could not fetch URL https://pypi.org/simple/laspy/: There was a problem confirming the ssl certificate: HTTPSConnectionPool(host='pypi.org', port=443): Max retries exceeded with url: /simple/laspy/ (Caused by SSLError(SSLCertVerificationError(1, '[SSL: CERTIFICATE_VERIFY_FAILED] certificate verify failed: unable to get local issuer certificate (_ssl.c:1129)'))) - skipping
ERROR: Could not find a version that satisfies the requirement laspy (from versions: none)
ERROR: No matching distribution found for laspy

処理が成功した際のQGISのverは3.34.5を使用しておりました。

（最新ver3.36や旧ver3.32では1度もlaspyのインストールが成功しませんでした...）

ありた様はどのverを使用していますでしょうか？

こういったエラー原因がご存知でしたら、ご教授頂けると助かります。

お手数をおかけしますが、ご確認お願い致します。

2024年5月24日金曜日 16:57:36 UTC+9 GTT:

[GTT]

ありた様

お世話になっております。

質問しておいてすいません、

エラーの原因を解消することが出来ました。

問題PCに設定されているセキュリティが問題でした。

プロキシを指定することで、laspyをinstallができました。（他のverでもinstall可能でした）

よろしくお願いいたします。

2024年5月27日月曜日 16:08:02 UTC+9 GTT: