ラベル 技術メモ の投稿を表示しています。 すべての投稿を表示
ラベル 技術メモ の投稿を表示しています。 すべての投稿を表示

2026年9月18日金曜日

技術メモ QGIS画面での分類結果データの書き出し方法

 Technical Note: Exporting Classification Results in QGIS


QGIS allows feature data to be displayed using various classification methods. I explored a way to export these classification results—represented by values ​​such as 1 through 10, depending on the number of classes—into a new data field. I created a convenient Python script that enables the export of this data regardless of the classification method or number of classes used in the QGIS display.


QGISでは地物データを多用な方法で分類表示できます。この分類表示結果を区分数に応じて1~10などの数値データで新たなフィールドに書き出す方法を検討しました。QGIS画面で表示したどのような分類方法、分類数でもデータ書き出しができる便利なPythonスクリプトを作成しました。

1 フィールド計算機を使った方法 その1

1-1 等間隔分類

サンプルデータ(0~90の範囲内にある小数点のある数値、レコード数3024、フィールドkousaに格納されている。)

フィールドkousaについて、シンボロジ→連続値による定義→フィールド入力→等間隔分類(分類数10)で画面上で分類しました。


画面上での分類結果

この画面上での分類結果を、数値の小さい方から1~10の数字を割り当てた区分データに変換し、それを新たなフィールドkousakubun3に出力する式(例)とフィールド計算機画面は次の通りです。

式(例)

CASE

    WHEN "kousa" IS NULL THEN NULL

    WHEN maximum("kousa") = minimum("kousa") THEN 1

    ELSE

        -- kousaの値を1.0〜10.999の範囲に線形変換し、floorで整数化(1〜10になる)

        floor(scale_linear("kousa", minimum("kousa"), maximum("kousa"), 1, 10.9999))END


フィールド計算機画面

これで無事等間隔分類を1~10の数値データに変換して、新たなフィールドに書き込むことができました。

1-2 等量分類

等量分類についての式は複雑になり過ぎて、ギブアップしました。

1-3 感想

等間隔分類はなんとかフィールド計算機で再現できましたが、等量分類はギブアップで、他の分類(丸め間隔、固定間隔、対数スケール、標準偏差、自然分類)は思考の意欲が生まれませんでした。

2 フィールド計算機を使った方法 その2

分類のアルゴリズムを再現するのではなく、QGIS操作で生じた分類結果そのものをフィールド計算機にとりこむ方法があります。

例えばシンボロジの分類結果が仮に、

1 :  2.14 ~  5.72

2 :  5.72 ~  8.31

3 :  8.31 ~ 11.45

...

10: 47.21 ~ 88.53

となった場合、フィールド計算機では、

CASE

WHEN "kousa" IS NULL THEN NULL

WHEN "kousa" <= 5.72 THEN 1

WHEN "kousa" <= 8.31 THEN 2

WHEN "kousa" <= 11.45 THEN 3

...

WHEN "kousa" <= 47.21 THEN 9

ELSE 10

END

とすれば、画面上の分類と同じ分類を属性として保存できます。

この方法なら、どのような分類でも、さらに一部を意図的に変更したものでも、対応できます。

ただし、分類の数値を手でひろい入力する必要があり、AI時代にふさわしい方法とは言えません。

3 Pythonスクリプトを使った方法

次のPythonスクリプトをプラグイン→Pythonコンソール→エディタの表示でひらくエディタ画面にはりつけて走らせれば、どのような分類でも1~10などの数値データを新たなフィールドのレコードに書き込むことができます。この方法も分類のアルゴリズムを再現するのではなく、QGIS操作で生じた分類結果そのものを新フィールドにとりこむ方法です。

このPythonスクリプトは、レイヤ名称aaa、対象フィールド名称bbb、新フィールド名称bbb_class、分類数10の場合で記述されていますが、それぞれ変更可能です。ChatGPT支援で作成しました。


from qgis.core import (
    QgsProject,
    QgsGraduatedSymbolRenderer,
    QgsField
)
from qgis.PyQt.QtCore import QVariant

# ============================================================
# 設定
# ============================================================

LAYER_NAME = "aaa"
VALUE_FIELD = "bbb"
CLASS_FIELD = "bbb_class"

# ============================================================
# レイヤ取得
# ============================================================

layers = QgsProject.instance().mapLayersByName(LAYER_NAME)

if not layers:
    raise Exception(f"レイヤ '{LAYER_NAME}' が見つかりません。")

layer = layers[0]

# bbb の存在確認
if layer.fields().indexFromName(VALUE_FIELD) == -1:
    raise Exception(f"フィールド '{VALUE_FIELD}' が見つかりません。")

# ============================================================
# 現在のシンボロジを取得
# ============================================================

renderer = layer.renderer()

if not isinstance(renderer, QgsGraduatedSymbolRenderer):
    raise Exception(
        "現在のシンボロジが「連続値による定義」ではありません。"
    )

# 分類に使用しているフィールドを確認
class_attribute = renderer.classAttribute()

if class_attribute != VALUE_FIELD:
    raise Exception(
        f"現在の分類フィールドは '{class_attribute}' です。"
        f"'{VALUE_FIELD}' ではありません。"
    )

# 現在画面で設定されている分類範囲を取得
ranges = list(renderer.ranges())

if len(ranges) != 10:
    raise Exception(
        f"現在の分類数は {len(ranges)} です。"
        "設定分類にしてから実行してください。"
    )

# ============================================================
# 分類範囲を確認表示
# ============================================================

print("現在の分類範囲")
print("------------------------------")

for i, r in enumerate(ranges, start=1):
    print(
        f"Class {i}: "
        f"{r.lowerValue()} ~ {r.upperValue()}"
    )

print("------------------------------")

# ============================================================
# bbb_class フィールドを作成
# ============================================================

class_idx = layer.fields().indexFromName(CLASS_FIELD)

if class_idx == -1:

    if not layer.isEditable():
        if not layer.startEditing():
            raise Exception(
                "レイヤを編集モードにできませんでした。"
            )

    success = layer.addAttribute(
        QgsField(CLASS_FIELD, QVariant.Int)
    )

    if not success:
        raise Exception(
            f"フィールド '{CLASS_FIELD}' を作成できませんでした。"
        )

    layer.updateFields()

    class_idx = layer.fields().indexFromName(CLASS_FIELD)

else:

    if not layer.isEditable():
        if not layer.startEditing():
            raise Exception(
                "レイヤを編集モードにできませんでした。"
            )

# ============================================================
# 各地物を現在のシンボロジと同じ区分に分類
# ============================================================

count = 0
null_count = 0
unclassified_count = 0

for feature in layer.getFeatures():

    value = feature[VALUE_FIELD]

    # NULLの場合
    if value is None or value == QVariant():
        layer.changeAttributeValue(
            feature.id(),
            class_idx,
            None
        )
        null_count += 1
        continue

    try:
        value = float(value)
    except (TypeError, ValueError):
        layer.changeAttributeValue(
            feature.id(),
            class_idx,
            None
        )
        unclassified_count += 1
        continue

    class_number = None

    # QGISの graduated renderer と同様に
    # 上から順番に該当する範囲を探す
    for i, r in enumerate(ranges, start=1):

        if (
            r.lowerValue() <= value
            and
            value <= r.upperValue()
        ):
            class_number = i
            break

    if class_number is not None:

        layer.changeAttributeValue(
            feature.id(),
            class_idx,
            class_number
        )

        count += 1

    else:

        layer.changeAttributeValue(
            feature.id(),
            class_idx,
            None
        )

        unclassified_count += 1

# ============================================================
# 結果
# ============================================================

print("")
print("処理完了")
print("==============================")
print(f"分類した地物数 : {count}")
print(f"NULL           : {null_count}")
print(f"分類不能       : {unclassified_count}")
print("==============================")
print("")
print(
    "属性テーブルの新フィールドに "
    "分類番号を書き込みました。"
)
print(
    "内容を確認して問題なければ、"
    "レイヤの編集内容を保存してください。"
)

4 感想


レイヤの属性テーブル 例

今回作成Pythonスクリプトでどのような分類、分類数でも、新しいフィールドのレコードに分類を指標する1~10など数値を書き込むことができるようになりました。画面の分類結果をデータとして残せるようになりました。画面で多種分類方法をシミュレーションして(分布を確認して)、適切なもののデータを残せるのでとても便利です。

2026年9月3日木曜日

技術メモ GeoPackageの中身を覗くPythonスクリプト

 Technical Note: Python Script to Inspect GeoPackage Contents


I created a Python script to inspect the contents of GeoPackage files (.gpkg) used in QGIS. When dealing with multiple GeoPackage files, this allows you to check the tables stored within them via the command line, without having to open each file individually in QGIS.


QGISで使うGeoPackageファイル(.gpkg)の中身を覗くPythonスクリプトを作成しました。GeoPackageファイルが多数ある場合、それらを一つ一つQGISで開かなくても、コマンドラインでGeoPackageファイルに格納されているテーブルを確認できます。

1 GeoPackageファイルに格納されている全てのテーブルを確認するPythonスクリプト

GeoPackageファイルに格納されている全てのテーブルの確認をコマンドラインで確認できるPythonスクリプトを作成しました。QGISで作成したテーブルだけでなく、SQLiteの管理用テーブルも含まれます。


GeoPackageファイルに格納されている全てのテーブルを確認している様子

2 GeoPackageファイルに格納されているQGIS関連テーブルを確認するPythonスクリプト

GeoPackageファイルに格納されているQGIS関連テーブルだけの確認をコマンドラインで確認できるPythonスクリプトを作成しました。SQLiteの管理用テーブルは含まれません。


GeoPackageファイルに格納されているQGIS関連テーブルを確認している様子

3 QGISPythonコンソールでGeoPackageファイルに格納されているQGIS関連テーブルと詳細情報を表示するPythonスクリプト

QGISPythonコンソールでGeoPackageファイルに格納されているQGIS関連テーブルと詳細情報を表示するPythonスクリプトを作成しました。詳細情報はレイヤ名・ジオメトリ型・CRS・フィールド名・型・地物数を一覧表として表示します。


QGISPythonコンソールでGeoPackageファイルの詳細情報を表示している様子

ただし、一覧表示でなければ、このPythonスクリプトで表示される情報はすべて通常のQGIS操作で確認できる情報です。

4 3のPythonスクリプトをコマンドラインで走らせるPythonスクリプト

3のPythonスクリプトをコマンドラインで走らせるPythonスクリプトを作成しました。


3のPythonスクリプトをコマンドラインで走らせるPythonスクリプトの様子

但し、現状では作業中のため、当該Pythonスクリプトと対象GeoPackageファイル(.gpkg)の双方のフルパス入力が必要なので、使い勝手が悪いです。

5 感想

GeoPackageファイルの内容はQGISにインストールすれば、通常操作で全て判ります。しかし、類似作業が増え、多数フォルダーに類似GeoPackageファイルが存在するようになった場合、それらを一つ一つQGISで確認することは非効率的です。多数GeoPackageファイルの確認に特に2のPythonスクリプトの存在意義が高まります。


2026年9月1日火曜日

技術メモ QGISのGeoPackage

 Technical Note: QGIS and GeoPackage


I am currently using QGIS to perform fabric analysis on 3D models of shell-layer profiles, and I frequently find myself copying the QGIS project file to a different folder for use. When working with the original QGIS file in a separate folder, storing the QGIS layers in a GeoPackage file is convenient because it means moving only two files in total—the main QGIS project file and the GeoPackage file.


貝層剥取断面3Dモデルのファブリック分析をQGISを使って行っていますが、同じQGISファイルを別フォルダーにコピーして使うことが多くなっています。別フォルダーで元QGISファイルを使う場合、QGISレイヤをGeoPackageファイルに格納すると、QGIS本体ファイルとGeoPackageファイルの合計2つのファイル移動だけで済みますから便利です。

1 GeoPackageを使わない場合の例

QGISで1つのshapeファイルの場合、shapeファイル本体(.shp)以外に.cpg,.dbf,.shxなどの関連ファイルが同時に必要になります。QGISファイルを別フォルダーにコピーして使うばあい、このshapeファイル1つだけのQGISでも、これら4つの本体・関連ファイルを同じフォルダーにコピーする必要があります。現実は多数レイヤになりますから、極めて多数のファイルを本体ファイルと一緒にコピーする必要があります。現実的対応ではありません。

2 GeoPackageを使う場合の例

多数レイヤ(ベクター、ラスター)を1つのGeoPackageに格納すれば、QGIS本体ファイルとGeoPackageファイルの2つだけを別フォルダーにコピーすれば、その別フォルダーで使えます。なお、QGIS本体のプロパティ→全般→保存パスを相対パスにしておく必要があります。(デフォルトで既にそうなっています。)

多数レイヤを1つのGeoPackageに格納する場合、プロセッシングツールボックス→レイヤをGeoPackage化を使えば、ベクターレイヤのシンボロジなど各種設定がそのまま保存される設定にデフォルトでなっています。


レイヤをGeoPackage化の画面(パラメーター)


レイヤをGeoPackage化の画面(ログ)


別フォルダーで元QGISファイルを開いて、レイヤの設定が保存されている様子

なお、ブラウザパネルのGeoPackageファイルにレイヤパネルのレイヤをdrag&dropで格納できますが、この場合シンボロジなどの各種設定は保存されません。

3 メモ

GeoPackageを自分と別人とのQGISデータのやりとりのツールと理解していました。しかしそれだけでなく、自分だけの作業でも別フォルダーで派生する類似関連作業のツールして気軽に使えることに気が付きました。QGISは特定フォルダーで集中して行うという昔の習慣が今頃になって改善しました。


2026年7月15日水曜日

技術メモ モニターに表示されている画像からのOCR Windows+Shift+T

 Tech Memo: OCR for On-Screen Images (Windows + Shift + T)


I discovered a Windows 11 keyboard shortcut (Windows + Shift + T) that allows for text extraction (OCR) from images displayed on the monitor—such as images on a webpage or text from unexpected error messages—and I have found it incredibly useful.


モニターに表示されている画像(web画面内の画像など)の中の文字テキストの抽出(OCR)や、突然表示されるエラー表示のテキスト抽出ができるWindows11ショートカットキー(Windows+Shift+T)があることを知り、大変便利に使いだしました。

1 Windows+Shift+T


使用例(Windows画面切手画像の文字のOCR)

結果

…………………

UNESCO


Patrimoine mondial 1997

…………………


使用例(エラー表示、Blender)

結果

…………………

レポート:エラー


×


Python: Traceback (most recent call last):

File "G:\パソコンとソフト\ChatGPT\20260625Blenderアドオン改良\UI版バージョン3開発\k9rtest22.blend\gpt_generated.py", line 54, in <module>

File "G:\パソコンとソフト\ChatGPT\20260625Blenderアドオン改良\UI版バージョン3開発\k9rtest22.blend\gpt_generated.py",line 7,in clear_scene

NameError: name 'bpy' is not defined

…………………

Windows+Shift+Tはとても便利です。PrintScreenキーからこの機能を利用することもできます。

参考サイト【Windows標準OCR】窓キー + Shift + T キー押してますか?

2 関連Windows11機能 Windows+V

クリップボード履歴活用ショートカットキーWindows+Vも便利です。過去のCtrl+Cでクリップボードに記憶した情報を最大25まで遡って閲覧して、再利用できます。なおテキスト・画像混載です。


Windows+V画面

参考サイト【Windows】窓キー+V、押してますか?

3 メモ

参考サイトは定期的メール「Qiita[キータ] <info@qiita.com>」の「先週いいねが多かった投稿ベスト20」からみつけたものです。「Qiita[キータ] <info@qiita.com>」というメールがどうして届くようになったのか記憶をたぐることができませんが、「先週いいねが多かった投稿ベスト20」は役立つ情報が多いです。

クリップボード履歴活用ショートカットキーWindows+Vと同じような機能をMS-DOS時代のVZエディターで利用していた記憶があります。大変便利だったのですが、Windows時代になりいつのまにか忘れていました。


2026年7月9日木曜日

技術メモ Blenderスプレッドシートパネルに遺物点群の属性データを表示する方法

 Technical Note: How to Display Artifact Point Cloud Attribute Data in the Blender Spreadsheet Panel


Using Geometry Nodes in Blender causes an issue where point cloud attributes are no longer displayed in the Object Data Properties. Therefore, I verified a method to display artifact attribute data from artifact point cloud datasets within the Blender Spreadsheet panel.


Blenderでgeometry nodesを使うと、オブジェクトデータプロパティに点群属性が表示されなくなる不都合が生じます。そこで、Blenderスプレッドシートパネルに遺物点群データの遺物属性データを表示する方法を確認しました。

1 Blenderスプレッドシートに遺物点群データの遺物属性データを表示する方法


Blenderスプレッドシートに遺物点群データの遺物属性データを表示する方法

ジオメトリノードで遺物点群が格納されたオブジェクト情報ノードにビューアーノードを接続すれば、スプレッドシートパネルに遺物属性データが表示されます。

この画面では約6万遺物、40属性データの一部が表示されています。

2 メモ

これまでオブジェクトデータプロパティ→属性で属性情報の有無をみてきました。しかし、geometry nodesをつかうとオブジェクトデータプロパティ→属性では属性情報が完全に表示されないことがあることを知りました。その不都合を回避するために、geometry nodesにおける属性情報確認はオブジェクト情報にビューアーを接続して、スプレッドシートパネルで確認することにしました。


2026年7月1日水曜日

技術メモ 遺物3D空間分析用geometry nodesの使いまわし方法

 Technical Note: Reusing Geometry Nodes for 3D Spatial Analysis of Artifacts

I have been making effective use of a Geometry Nodes setup to analyze 3D spatial data for approximately 56,000 artifacts from the shell midden on the northern slope of the Ariyoshi-Kitakaizuka Shell Mound. To enable the reuse of this setup across different Blender files, I have converted it into an asset.


有吉北貝塚北斜面貝層の遺物約56000件データを3D空間で分析するgeometry nodesを便利に使っています。このgeometry nodesを別Blenderファイルでも使えるようにアセットにして使いまわし出来るようにしました。

1 geometry nodesを使っている様子


geometry nodesを使ている様子

2 遺物3D空間分析用geometry nodesのイメージ


遺物3D空間分析用geometry nodesのイメージ

このgeometry nodesには次の機能が備わっています。

1)3D座標付遺物データの選択

このgeometry nodesでは有吉北貝塚北斜面貝層の全42属性、レコード数55892件のデータを選択してあります。

2)表示する属性の選択

属性bunruicodeを選択すれば、それは遺物台帳に紐づくデータで全遺物が1(土器)、2(石器)、…のように区分されています。

属性zinkotsuを選択すれば、それは発掘調査報告書に紐づくデータで、該当データだけ区分が1になっています。

3)表示する属性の区分の選択

属性bunruicodeの区分1は土器を、区分5は骨歯を選択します。

属性zinkotsuの区分1は発掘調査報告書掲載人骨データを選択します。

4)表現メッシュ形状の選択

vert(点[メッシュ])の他、球、立方体など好みに合わせて表現メッシュ形状を選択できます。球や立方体は大きさや色なども設定できます。

3 geometry nodesの使い勝手

geometry nodesでノード(箱)の中の区分数値を変化させるだけで、表現内容が一瞬のうちに変化しますので、とても使い勝手がよいです。

しかし、このgeometry nodesのノード(箱)を一つ一つ組み立て、ラインで結ぶのは根気のいる作業です。

4 geometry nodesの使いまわし

別のBlenderファイルでも遺物3D空間分析用geometry nodesを使いまわす方法は幾つかありますが、これまではアペンドを使ってきました。しかしアペンドは遺物3D空間分析用geometry nodesが存在するBlenderファイルから持ってくる方法なので、このgeometry nodesが存在しているBlenderファイルを探さなければなりません。意外と手間取ります。

今回geometry nodesのアセット登録してみると、とても効率的にgeometry nodesを使いまわせることを確認しました。


アセットによるgeometry nodesの使いまわし方法

アセットによるgeometry nodesの使いまし方法の概要は以下のとおりです。

1)geometry nodesのグループ化

2)グループを「アセットとしてマーク」→これでアセットが使えるようになります。(※)

3)別Blenderファイルでplane追加

4)planeに新期geometry nodes作成して、そこにアセットをドロップして、グループとグループ出力をラインで結ぶ

これで遺物3D空間分析ができるようになります。(なお、アセットをドロップすると、元遺物データとvertが自動生成します。)

※ 「アセットとしてマーク」を実行するファイルを事前に設定するアセット専用フォルダーに置いておく必要があります。

2026年6月30日火曜日

技術メモ BlenderPythonスクリプト資産の効率活用策

 Technical Note: Strategies for Efficiently Utilizing Blender Python Script Assets


As part of an effort to efficiently utilize the library of Blender Python scripts I have built up over time, I have adopted the practice of making appropriate use of the File Browser panel within the Blender interface. This has facilitated the efficient use of existing scripts.


これまでの活動で貯めてきたBlenderPythonスクリプト資産の効率的活用策の一環としてBlender画面でファイルブラウザーパネルを適宜使うようにしました。これにより既存スクリプトの効率的活用が進んでいます。

1 Blender画面におけるファイルブラウザーパネルの設置とBlenderPythonスクリプト資産の表示


BlenderPythonスクリプト資産を表示したファイルブラウザーパネル設置Blender画面

パソコンに整理されているBlenderPythonフォルダーをブラウザパネルに表示して、BlenderPythonスクリプトを選びやすくしています。

2 Blender画面テキストエディターにおけるBlenderPythonスクリプト説明文書の表示


BlenderPythonスクリプト説明文書を表示しているBlender画面

BlenderテキストエディターはBlenderPythonスクリプトだけではなく、説明文書も表示できますから、BlenderPythonスクリプト資産活用の際にスクリプト内容についてより深く確認できます。

3 テキストエディターのBlenderPythonスクリプトを走らせた様子


テキストエディターでBlenderPythonスクリプトを走らせた画面

パソコンフォルダーに格納してあるBlenderPythonスクリプトをファイルブラウザー経由で取り出しテキストエディターで走らせることによって、エクスプローラー経由よりも効率的に作業が進みます。

4 メモ

これまでWindowsのエクスプローラーを使って既存資産を取り出し、Blenderに貼りつけていましたが、それより、ファイルブラウザーを使う方が効率的であることを実感しています。効率化実感の一因はウィンドウ数減少にあります。


2026年6月22日月曜日

技術メモ pngファイルからヒートマップを作成するPythonスクリプトの改良

 Technical Notes: Improvements to a Python Script for Creating Heatmaps from PNG Files


I've been using a very convenient Python script to create heatmaps from dot distribution images (transparent background PNG files). However, I encountered some issues, so I've improved it. The resulting Python script is now more robust.


ドット分布画像(背景透明pngファイル)からヒートマップを作成するPythonスクリプトをとても便利に使っています。しかし、不都合が生じたので改良しました。より堅牢なPythonスクリプトとなりました。

この記事は2026.06.03記事「技術メモ pngファイルから直接ヒートマップを作成するPythonスクリプト」の続きになります。

1 遭遇した不都合


遭遇した不都合

これまで順調にpngファイルからPythonスクリプトで直接ヒートマップを作成してきました。しかし、上記のような不都合に遭遇しました。

調べたところ、画像の透明部分にRGBの色が残っている部分が存在しているために生じた現象であることがわかりました。Pythonスクリプトは透明か否かではなく、色がある部分=点群データとして扱っているので生じた不都合です。

2 不都合の解決策

この不都合の解決策として、次の2を執りました。

1 画像透明部から色を完全に除去する新たなPythonスクリプト作成。

2 色のある透明部を間違って拾わないヒートマップ作成Pythonスクリプトの作成。

それぞれのPythonスクリプトを作成しましたが、2の新たな改良ヒートマップ作成Pythonスクリプトを使う方が作業の手間がはぶけます。ヒートマップ作成Pythonスクリプトがより堅牢となりました。

3 改良Pythonスクリプト (透明部に色が残っていても正常に作動する改良版)


from pathlib import Path

import cv2
import numpy as np
from PIL import Image
from matplotlib import cm

# =========================
# 設定
# =========================

input_path = r"G:/test/aaa.png"
output_path = r"G:/test/aaa_heatmap.png"

# 使用する単色グラデーション
# 例: "Blues", "Reds", "Greens", "Purples", "Oranges", "Greys"
color_map_name = "Reds"

# グラデーションの区分数
# 例: 5, 8, 10, 16, 32
gradient_steps = 10

# 密度計算のぼかし半径
# 大きいほど広くなめらかな密集度になる
blur_radius = 20

# 背景を透明にするか
transparent_background = False

# 点群抽出に使うアルファ閾値
# alpha > alpha_threshold のピクセルだけを点として扱う
# 0: わずかでも不透明なら点
# 10〜30: ほぼ透明なノイズを除外
alpha_threshold = 20

# ヒートマップの透明出力時、薄い密度部分を透明にする閾値
# 通常は 0.0 のままでOK
heat_alpha_threshold = 0.0

# =========================
# 処理
# =========================

input_file = Path(input_path)
if not input_file.exists():
    raise FileNotFoundError(f"画像が読み込めません: {input_path}")

# RGBAで読み込み、RGB値には依存しない
img = Image.open(input_file).convert("RGBA")
arr = np.array(img)

# アルファチャンネルだけで点群を抽出
alpha = arr[:, :, 3]
dot_mask = alpha > alpha_threshold

point_pixels = int(np.count_nonzero(dot_mask))
if point_pixels == 0:
    raise ValueError(
        "アルファ値から点を検出できませんでした。"
        f" alpha_threshold={alpha_threshold} を下げてください。"
    )

print(f"画像サイズ: {img.size[0]} x {img.size[1]}")
print(f"alpha min/max: {int(alpha.min())} / {int(alpha.max())}")
print(f"点として使うピクセル数: {point_pixels}")

# 密度画像を作成
density = dot_mask.astype(np.float32)

# ガウシアンぼかしで密集度を計算
density = cv2.GaussianBlur(
    density,
    ksize=(0, 0),
    sigmaX=blur_radius,
    sigmaY=blur_radius,
)

# 0〜1に正規化
max_density = float(density.max())
if max_density > 0:
    density = density / max_density

# グラデーションを段階化
density_step = np.floor(density * gradient_steps) / gradient_steps
density_step = np.clip(density_step, 0, 1)

# カラーマップ適用
cmap = cm.get_cmap(color_map_name)
heat_rgba = cmap(density_step)

# 0〜255へ変換
heat_img = (heat_rgba[:, :, :3] * 255).astype(np.uint8)

# 密度ゼロ部分の処理
if transparent_background:
    out_alpha = (density > heat_alpha_threshold).astype(np.uint8) * 255
    output = np.dstack([heat_img, out_alpha])
    Image.fromarray(output, mode="RGBA").save(output_path)
else:
    # 背景は白
    heat_img[density <= heat_alpha_threshold] = [255, 255, 255]
    Image.fromarray(heat_img, mode="RGB").save(output_path)

print(f"ヒートマップを書き出しました: {output_path}")

2026年6月3日水曜日

技術メモ pngファイルから直接ヒートマップを作成するPythonスクリプト

 Technical Memo: Python Script to Create Heatmaps Directly from PNG Files


I created a Python script that directly outputs heatmaps from dot distribution images (PNG files). Previously, I used QGIS functions to create heatmaps, but this significantly improves work efficiency.


ドット分布画像(pngファイル)から直接ヒートマップをアウトプットするPythonスクリプトを作成しました。これまではQGIS機能を利用してヒートマップを作成していたのですが、それと較べると著しく作業効率が向上しました。

1 Pythonスクリプトによるヒートマップ作成


利用画像(ドット分布画像、背景が透明なpng画像)

ドット分布画像を用意します。背景が白地の場合は背景が透明なpngファイルにします。

ドットの色の面積をpixel単位に計測して、それによりヒートマップを作成するので、ドットが重なってつぶれていても、データは近似的に利用できます。


ヒートマップ

単色グラデーションによるヒートマップを作成します。

ヒートマップ作成は次の項目を変更して望みのイメージに合うまで試行できます。

色、グラデーションの区分、ぼかし半径


ヒートマップとドット分布図のオーバーレイ画像

2 使い勝手

これまではQGIS機能を利用してヒートマップを作成していたのですが、Pythonを使うと著しく作業効率が向上しました。

ドット分布図のドットが重なってつぶれているような場合でも、データとして近似的に利用できるので、ヒートマップの趣旨(分布概要をざっくり把握する)とも合います。

3 Pythonスクリプト


import cv2
import numpy as np
from PIL import Image
from matplotlib import cm

# =========================
# 設定
# =========================

# ドット分布図のパスとアウトプットのパス→変更してください
input_path = r"G:/test/aaa.png"
output_path = r"G:/test/aaa_heatmap.png"

# 使用する単色グラデーション
# 例: "Blues", "Reds", "Greens", "Purples", "Oranges", "Greys"
color_map_name = "Purples"

# グラデーションの区分数
# 例: 5, 8, 10, 16, 32
gradient_steps = 10

# 密度計算のぼかし半径
# 大きいほど広くなめらかな密集度になる
blur_radius = 50

# 背景を透明にするか
transparent_background = False

# =========================
# 処理
# =========================

# 画像読み込み
img = cv2.imread(input_path, cv2.IMREAD_COLOR)
if img is None:
    raise FileNotFoundError(f"画像が読み込めません: {input_path}")

# RGBへ変換
rgb = cv2.cvtColor(img, cv2.COLOR_BGR2RGB)

# 黒背景以外の点を抽出
# 青い点など、黒でないピクセルを遺物ドットとみなす
gray = cv2.cvtColor(rgb, cv2.COLOR_RGB2GRAY)
dot_mask = gray > 1

# 密度画像を作成
density = dot_mask.astype(np.float32)

# ガウシアンぼかしで密集度を計算
density = cv2.GaussianBlur(
    density,
    ksize=(0, 0),
    sigmaX=blur_radius,
    sigmaY=blur_radius
)

# 0〜1に正規化
if density.max() > 0:
    density = density / density.max()

# グラデーションを段階化
density_step = np.floor(density * gradient_steps) / gradient_steps
density_step = np.clip(density_step, 0, 1)

# カラーマップ適用
cmap = cm.get_cmap(color_map_name)
heat_rgba = cmap(density_step)

# 0〜255へ変換
heat_img = (heat_rgba[:, :, :3] * 255).astype(np.uint8)

# 密度ゼロ部分の処理
if transparent_background:
    alpha = (density > 0).astype(np.uint8) * 255
    output = np.dstack([heat_img, alpha])
    Image.fromarray(output, mode="RGBA").save(output_path)
else:
    # 背景は白
    heat_img[density <= 0] = [255, 255, 255]
    Image.fromarray(heat_img, mode="RGB").save(output_path)

print(f"ヒートマップを書き出しました: {output_path}")

2026年5月28日木曜日

技術メモ 型式別土器3D分布モデルの改善

 Technical Memo: Improvement of the 3D distribution model of Pottery by Type


The 3D distribution model of pottery by type was created in 2023 when the technology was still immature, resulting in data that was difficult to use. Therefore, I spent a day refining the old data to make it usable with Blender's geometry nodes. Simultaneously, I created a Blender Python script to project the point cloud data selected in geometry nodes onto a neighboring cross-section.


型式別土器3D分布モデル作成は技術が未熟であった2023年作業であるため、使い勝手の悪いデータとなっていました。そこで古いデータを1日かかりでBlenderのgeometry nodesで使えるデータに仕立て上げました。同時に、geometry nodesで選択した点群データを近傍断面図に投影するBlenderPythonスクリプトを作成しました。

1 型式別土器3D分布モデル(点群データ)の改善

旧データは土器片1件1件が独立したCUBEオブジェクトとなっています。型式別に分けて整理していますが、1件1件のデータには土器型式情報が含まれていません。

CUBEオブジェクトを球にしたり、大きさを変更したり、型式別に色を付けるなどの操作ができないことはないのですが、操作自由度がとても低くなっています。

今回、ID、XYZ座標のほか、土器型式分類の属性を加えたcsvファイルを作成し、属性を含めてBlenderにとりこみました。

csvファイルの項目:ID、X、Y、Z、bunrui(型式区分 1~9)

土器型式区分は次の通りです。

1 早期土器

2 前期土器

3 阿玉台式土器

4 中峠式土器

5 加曽利EⅠ式土器

6 加曽利EⅡ式古段階土器

7 加曽利EⅡ式中段階土器

8 加曽利EⅡ式中~新段階土器

9 加曽利EⅡ式新段階土器

2 属性(型式区分)を保持してBlenderにインポート

属性(型式区分)を保持してcsvファイルをBlenderにメッシュオブジェクト(点群)としてインポートしました。インポートはPythonで行いました。このインポートした点群データはgeometry nodesで、土器型式区分別に表現したり、表現するオブジェクトの形・大きさ・色などを自由に変更できます。


geometry nodesの様子


Blender画面

3 遺物投影断面図作成Pythonスクリプトの作成

既にBlender3Dビューポートに配置されている断面図を軸に前後0.5mのボックスを設定し、そのボックス内の遺物を抽出し、抽出遺物と断面図をオーバーレイして画像としてエクスポートするBlenderPythonスクリプトを作成しました。

対象となる遺物はgeometry nodesで抽出した遺物も対象にすることができます。

従って、このBlenderPythonを使うことにより、次のような遺物投影断面図を作成することができます。

型式別土器3D分布モデルからgeometry nodesで「7 加曽利EⅡ式中段階土器」だけの3D分布が表現されている状況があるとします。その「7 加曽利EⅡ式中段階土器」を特定断面に投影した投影断面図を作成することができます。投影する幅や投影断面図の表現(ドットの大きさや色)は自由に設定できます。

このBlenderPythonスクリプトを使うことにより、土器型式と貝層分層の関係把握を前進させることができると期待します。


第6断面に全遺物を投影させた様子(点メッシュを投影)


同部分拡大図


第6断面に「7 加曽利EⅡ式中段階土器」を投影させた様子(球を投影)


同部分拡大図


2026年5月20日水曜日

技術メモ QGISの疑似座標GIS解析を使った等高線からの標高点群生成

 Technical Memo: Generating Elevation Points from Contour Lines Using QGIS's Pseudo-Coordinate GIS Analysis


I have noted down a technique for generating elevation points from contour lines drawn on a plan view of the shell bed on the north slope of the Ariyoshi Kita Shell Mound. This technique utilizes QGIS's pseudo-coordinate GIS analysis function.


有吉北貝塚北斜面貝層の平面図に記入されている等高線から標高点群を生成する技術をメモしました。QGISの疑似座標GIS解析機能を利用したものです。

1 有吉北貝塚北斜面貝層平面図のQGISジオリファレンス


有吉北貝塚北斜面貝層平面図(QGISジオリファレンス用メモ記入)

1-1 QGISで使うEPSG

QGISのEPSG6677(日本測地系2011(JGD2011)[平面直角座標系]のⅨ(9)系(東京都 福島県 栃木県 茨城県 埼玉県 千葉県 群馬県 神奈川県)を利用しました。

1-2 座標変換

平面図にはもともとBlender投影用の座標が設定されています。Blender座標体系とQGIS座標体系は異なるので、またQGIS仕様に対応するため、今回次の座標変換等の対応と操作をしました。

・Blender用座標原点(0,0)とQGIS(EPSG6677)原点(0、0)を一致させる。

・Blender用座標X軸の値を正負逆転させてから1000倍して、QGIS座標Y軸に変換する。

・Blender用座標Y軸の値を1000倍して、QGIS座標X軸に変換する。

値を1000倍したのは、Blender用座標値そのままでは値が小さすぎてQGISが動かないためです。

1-3 ジオリファレンス

今回は5座標の値(-2000,0)( 0,0)( 40000,0)( -20000,-20000)( 40000,-20000)でジオリファレンスしました。


有吉北貝塚北斜面貝層平面図のQGISジオリファレンス結果(背景マップは作業と無関係です)

2 等高線トレース

2-1 新規シェープファイル作成

新規シェープファイルレイヤをファイル名とともに、ジオメトリ型:ラインストリング、追加次元:Z値、新規属性名称:elev、新規属性型:倍精度不動小数点型で作成します。

2-2 トレース

編集モードに入り、等高線をトレースします。1本1本についてトレース完了後、ID(通し番号)とelev(標高値)を記入します。


等高線トレース結果(背景マップは作業と無関係です)

3 標高点群生成とcsvファイル出力

3-1 標高点群生成

等高線レイヤを選択の上、プロセッシングツールボックスから「ジオメトリに沿った等間隔点群」を起動します。

ここでは距離を1000m(間隔)にしました。(平面図上では1m間隔に該当します。)

実行すると標高点群が生成します。


生成した標高点群(背景マップは作業と無関係です)

3-2 XYZフィールドの生成

生成した標高点群には標高はelevフィールドで表現されていますが、XY座標フィールドがありません。そこで、フィールド計算機を起動して、XYZフィールドを作成します。

Xフィールドは式に$Xを、Yフィールドは式に$Yを、Zフィールドは式に"elev"を入力してそれぞれのフィールドを生成します。


XYZフィールドが作成された属性テーブルの様子

3-3 csvファイル出力

XYZフィールドが作成された標高点群レイヤをcsvファイルで出力します。

4 標高点群からBlender用地形3Dモデル作成

4-1 csvファイルの調整

csvファイルをBlender用座標に戻します。

・Y座標を正負逆転させてから1/1000倍してX座標とする。。

・X座標を1/1000倍してY座標とする。

4-2 Blender用地形3Dモデル作成

今回は等高線から作成した標高点群を既存標高点群とマージして、それからBlenderPythonで地形3Dモデルを作成しました。

過去には、疑似座標GIS解析でトレースした等高線や標高点群からQGIS機能(v.surf.rstツール)で地形3Dモデルを作成したこともあります。

等高線や標高点群から地形3Dモデルを作成する方法はさらにいろいろあるので、それぞれを試して、使い勝手の良いものを見つけることにします。


2026年4月27日月曜日

技術メモ 貝層剥取断面と二枚貝(ハマグリ)の交線描画原理

 Technical Memo: Principle of Drawing Intersection Lines Between Shell Layer Cross-Sections and Bivalve Shells (Clams)


To measure the inclination between shell layer cross-sections and bivalve shells (clams), it is necessary to draw the intersection lines between the cross-section and the shell. This memo illustrates the principle of this drawing method.


貝層剥取断面と二枚貝(ハマグリ)の傾斜を計測するためには断面と貝殻の交線を描く必要があります。その描画方法の原理を図解でメモしました。

1 検討例


検討例

上凸で傾いたハマグリと下凸で傾いたハマグリの2例について、貝層剥取断面(鉛直平面)とハマグリの交線描画法原理を以下にメモします。

2 上凸で傾いたハマグリの貝殻傾斜線描画原理


上凸で傾いたハマグリ例

貝殻縁に密着する仮想平面をここでは殻縁平面と呼ぶことにします。殻縁平面と鉛直平面の交線が貝殻傾斜線となります。貝殻傾斜線と水平面との成す角度が貝殻傾斜となります。貝殻傾斜線描画後の傾斜角度計測はアプリ(自作Pythonスクリプト)が自動で行います。


実際例

実際の貝殻傾斜線描画は上記原理モデルを念頭に直観的に傾斜線を想定して、その端点2ヶ所をクリックして行います。

なお、傾きの状況が判りずらい場合は3Dモデルを観察して確かめてから描画します。

3 下凸で傾いたハマグリの貝殻傾斜線描画原理


下凸で傾いたハマグリ例


実際例

4 メモ

ここで求める貝殻傾斜はあくまでも剥取断面に対する傾斜であり、地山地形最大傾斜に対応する断面に対する傾斜ではないので、利用上注意が必要です。

また描画原理に基づくとはいえ、実際の作業は直観的な作業になりますから、精度の点では期待度の低いデータとなります。しかし、このようなデータは他にないので、貴重なデータとなります。


2026年4月26日日曜日

技術メモ 貝層剥取断面画像に見える地物大きさの視認性を向上させる画像処理

 Technical Memo: Image Processing to Improve the Visibility of Surface Object Size in Shell Layer Cross-Sectional Images


Through trial and error, I have acquired an image processing technique that significantly improves the visibility of surface objects such as seashells in shell layer cross-section images. This allows for intuitive understanding that, within the same shell layer, surface objects are larger at the lower part of the slope and smaller at the upper part.


貝層剥取断面画像に見える貝殻など地物大きさの視認性を大幅に向上させる画像処理技術を試行錯誤のなかで獲得しました。同一貝層で、斜面下部では地物大きさが大きく、斜面上部では小さいことを直観的によく理解できます。

1 貝層剥取断面画像


貝層剥取断面画像 (3Dモデルのオルソ投影画像)

貝層分層線は有吉北貝塚北斜面貝層セクション図分層線を画像に合わせて投影トレースしたものです。

2 地物大きさの視認性向上画像


地物大きさの視認性向上画像

貝層剥取断面画像をPhotoshopで「輪郭検出」したものです。

3 参考 セクション図における貝層分層


参考 セクション図における貝層分層

有吉北貝塚北斜面貝層セクション図の貝層分層記号です。分層記号の記述(説明)は掲載されていません。

4 メモ

地物大きさの視認性向上画像のC(2)層を見ると、斜面下部では地物大きさが大きく、斜面上部では小さいことを直観的によく理解できます。

C層(上)、C(2)層は二枚貝(ハマグリ)をメインとする貝層、C層(下)、C(3)層はイボキサゴをメインとする貝層です。


2026年3月2日月曜日

技術メモ 展示物の3Dモデル作成用撮影について

 Technical Note: Photographing Exhibits for 3D Modeling


When I photographed a double-mouthed pottery vessel on display in the special exhibition "The Mysterious Late Jomon Period" and created a 3D model, the 3D model completely failed. However, after reworking it with just one additional photograph, I was able to create a satisfactory 3D model. I investigated the reasons for this and made notes.


特別展「謎多き縄文晩期」で展示されている双口土器を周回撮影し、3Dモデルを作成したところ、3Dモデルが完全に破綻しました。しかし、たった1枚の写真を追加して再作業すると、立派な3Dモデルができましたので、その理由を検討し、メモしました。

1 3Dモデル作成 ケース1

千葉市埋蔵文化財調査センター特別展「謎多き縄文晩期」で展示されている弥生前期双口土器の3Dモデルを作成しました。170枚の写真を多段(6段)周回撮影しました。


多段周回撮影の様子

撮影は双口土器が画面の中央にくるように撮影しました。


双口土器撮影写真例


撮影写真

フォトグラメトリソフト3DF Zephyr Liteで3Dモデルを作成したところ、完全に破綻した3Dモデルとなりました。


完全に破綻した3Dモデル

3DF Zephyr Liteがはねた写真は9枚です。

2 3Dモデル作成 ケース2

ケース1で使った170枚写真の先頭に次の展示台全体が写る写真1枚を加えました。


ケース1の170枚に加えた1枚

3DF Zephyr Liteで同じ作業をしたところ、次の完全に正常な3Dモデルができました。


正常に造形された3Dモデル

3DF Zephyr Liteがはねた写真は1枚です。

3 写真1枚で3Dモデル造形が劇的に改善した理由

ケース1に写真1枚を加えたケース2で、3Dモデル造形が劇的に改善した理由について、ChatGPTの支援をうけて検討し、次のように想定しました。

3-1 ケース1で3Dモデルが破綻した理由

・土器が黒く微細テクスチャに乏しく、土器以外がほぼ白く無地で安定した特徴点(対応点)が得られない。

・そのため誤対応が混ざり、誤った姿勢推定となり、メッシュ化で破局して、紙が剥がれたようになった。

3-2 1枚の写真で造形が劇的に改善した理由

・加えた1枚の写真には展示台、壁の継ぎ目、説明板、他遺物、ラベルなど形がはっきりした特徴点が大量に増え、カメラ姿勢が一気に安定した。これは結果として「はねられる写真」が 9→1 に減ったこととも整合する。

・加えた1枚の写真を先頭に入れることによって、スコアの高い特徴点が採用され、初期解が最適化した。3DF Zephyr LiteのようなSfMでは、最初の数枚で作る骨格が悪いと、その後ずっと悪い解を引きずる。

・結果として1枚で劇的に改善された。

・加えた1枚の写真は構図が広いので、画面全域に直線や構造物が入りやすく、歪み推定の手掛かりが増えるため、内部パラメータが安定→姿勢推定も安定、につながった。

4 メモ

展示物3Dモデル作成では150枚程度撮影するとそれなりの品質のものができることを経験的につかんできています。しかし、今回のように真っ黒で表面テクスチャに特徴がなく、かつ周辺に別の地物がないと、3Dモデルが破綻することを改めて理解しました。

今後、展示物3Dモデル作成用撮影では構図の広い写真を数枚以上必ず意識的に追加することとします。

この双口土器は次の記事で説明してます。

2026.02.17記事「弥生時代前期 双口土器(富岡市上高田社宮子原遺跡)観察記録3Dモデル



2026年1月31日土曜日

技術メモ GigaMesh Software FrameworkにおけるNon-Photorealistic Renderingの表現工夫

 Technical Note: Techniques for Using Non-Photorealistic Rendering in the GigaMesh Software Framework


When creating 3D models of Jomon pottery, I usually use the GigaMesh Software Framework to flatten the image so that the patterns can be observed at a glance. When using non-photorealistic rendering for these flattened images, I found that adding coloring made them more visible.


縄文土器3Dモデルを作成する時は、ほとんどの場合GigaMesh Software Frameworkで平面展開して、文様を一覧的に観察できるようにしています。この展開画像としてNon-Photorealistic Renderingを使う場合、彩色することでより見やすくできることがわかりました。

この記事は2026.01.30記事「浅鉢(成田市大原野貝塚)観察記録3Dモデル」の続きです。

1 GigaMesh Software FrameworkによるNon-Photorealistic Rendering表現


GigaMesh Software Frameworkへの土器3Dモデルのインポート


GigaMesh Software Frameworkでの展開(Solid)


GigaMesh Software Frameworkでの展開(Non-Photorealistic Rendering)

線分成分だけを彩色(ここでは赤色)し、他の成分は除去しました。

2 表現工夫


Solid画像


Non-Photorealistic Rendering画像


Solid画像とNon-Photorealistic Rendering画像の乗算ミックス画像

Solid画像とNon-Photorealistic Rendering画像の乗算ミックス画像を作成すると、文様の様子が一段と強調され、判りやすくなります。

3 メモ

Non-Photorealistic Rendering画像は世間一般では漫画に使ったり、キャラクター表現に使ったりなど、正確性・写実性から退化させるための用途で使われています。しかし、GigaMesh Software Frameworkでは正反対に、粘土板の象形文字とか土器の文様分析など表面の細かいディテールを可視化する手法として用いられています。今後縄文土器文様や縄文施文のディテール分析に有効なツールとして使っていくことにします。