はじめに

広義の「動物」（哺乳類だけでなく、ここでは鳥類も昆虫類も爬虫類などもふくみます）は、太古の昔より自ら移動することで栄養を得て、子孫を残してきました。そして研究者たちはこの「動物」の「移動」の不思議…つまり「なんでアンタそこにいるのよ？」の解明にのために様々な工夫や調査手法を開発してきました。
その古典的手法のひとつがVHFテレメトリー法という調査で、VHF電波発信機を動物に装着し、個体の足取りや環境や生活史とのかかわりを考察してきました。（詳細は佐伯・早稲田(2006)などをお読み下さい）。

via. wikipedia

さて、この手法の特徴は簡易三角測量によって対象個体の位置推定を行うことでした。つまり、電波感度の良い測位点(A-C)から指向性の高いアンテナで個体位置(🐻)の方角を推定し、コンパスで方位を測定し、可能な限りすばやく別の地点から同じ作業を行い、地形図に手書きで線を引くことでした。様々な手法があるものの、可能な限りこの三角形が小さいほうがよいとされ、その重心点が動物の推定位置とみなしてデータ収集を行っています。

そしてさらに大変なところは、研究者たちはこの作業を何回も、何日も、そして場合によっては何年も続けなくてはいけないということです。

でも、令和のヤングたちはこう思うでしょう
👨‍🎓「なんでGPSを使わないんすかー？」
そうですよね、GPSがこの分野でも心おきなく利用できるようになったのは2010年代前後であり、重量・コスト・精度・頻度などによるメリット・デメリットなどの制約から、いまでもGPSとVHFの手法は共存している状況です（たとえばDore et al., 2020)。

わたしも学生時代からVHFテレメトリを用いて卒論も書いたり、ヤマネコや何十kmも長距離移動を行うシカやクマの位置をさぐる調査プロジェクトに参加したりしていました。そんなバックグラウンドもあって「なんでアンタそこにいるのよ？」の観点や好奇心は今も変わらず、これをもとに様々なお仕事や社会貢献のお手伝いができているのは本当にありがたいことです。

↑ヒグマ捕獲に携わる若きわたし（ケツ向けてますが）

とはいえ、VHFテレメトリのデータ処理を何千回と手作業で行うことは苦行でした。上記の各種論文でも触れられているのですが、VHF方式は重量や価格といったハードウェアの制約が小さい分、データ取得や処理に関わる人的コストが莫大に大きいという課題があります。その後GISの発達などにより様々な前処理ツールも出現しましたが、未だに方向性が定まっていない状況にあります。

そこでこの記事では、VHFテレメトリ法を利用し、今もなおフィールドでもがき苦しむ研究者や学生さんに向けてQGISでのVHFデータ前処理について調査を行い、その処理方法を提案したいと思います。

データ準備

今回、この検証を行うために、日本国内のVHFテレメトリー法による生データを探しましたがたどり着けませんでした。そこで、世界有数の動物測位データポータルMovebankから拝借しようと思ったのですが、ここでも日本国内の連続的な良いデータはヒットしませんでした(ほぼ渡り鳥とかコウモリとか…)。

そこで、大学一年のとき無人島で海水茹でパスタを一緒に食って🤮した渡辺くんさんのご尽力によりBiP-EarthにてCC-BY公開された、山口県農林総合技術センター松本哲朗氏によるアライグマのGPSトラッキングデータから、VHFテレメトリダミーデータをスクリプトで作成しました(a.k.a #餅から米)

GPS装着個体と追跡調査風景（渡辺氏提供）

VHFテレメトリー法で得られるデータは基本、いつ、どのセットで、どこで、どの方向に、どの個体がいた？となり、今回はこのデータ構造をベースとします。こんな感じですね。

カラム名	内容	例
name	個体の識別子	PL018H
set_id	観測セット番号	3
datetime	測位日時	2021-01-23 10:25:00
lat	観測地点の緯度（WGS84）	34.5647877
lon	観測地点の経度（WGS84）	131.6648722
azi_deg	測定した方位角（度）	275.6

なお、前述のようにGISもなにもない時代は私を含め多くのフィールドワーカーは全部これを地形図に手で線を引いていました。ﾋ-

さて、今回用いるダミーデータは下記になります。これをvhf_result.csvなどとして任意の場所に保存してください。アライグマ２個体の測位点10点について、道路からVHFテレメトリーで各々３回、合計30回測位した場合のダミーデータです。
時間はランダムですし、道路上での測位点も地形などを考慮していない実情に沿わないデータになっていますが、ご容赦ください。また、今回は磁気偏角は考慮していませんのでご注意を。

name,set_id,datetime,lat,lon,azi_deg
PL018H,1,2021-01-23 09:11:00,34.57453220894295,131.66278412612684,177.3
PL018H,1,2021-01-23 09:23:00,34.55568911004033,131.6519302429112,61.5
PL018H,1,2021-01-23 09:40:00,34.55188604321555,131.6806039156737,305.7
PL018H,2,2021-01-23 10:28:00,34.5471253091056,131.64277397144317,21.3
PL018H,2,2021-01-23 10:42:00,34.56299806374202,131.66498515236833,249.6
PL018H,2,2021-01-23 10:59:00,34.56313212021094,131.63515247910593,125.1
PL017H,1,2021-01-23 11:34:00,34.56222225723217,131.66486620903012,23.5
PL017H,1,2021-01-23 11:45:00,34.57305130529924,131.68034459642175,276.8
PL017H,1,2021-01-23 12:01:00,34.58466598254374,131.66220060185967,149.3
PL017H,2,2021-01-23 13:19:00,34.55157225697689,131.6557436056433,91.5
PL017H,2,2021-01-23 13:38:00,34.5379835718385,131.6815466983977,339.6
PL017H,2,2021-01-23 13:48:00,34.55900763780229,131.68401127645044,228.7
PL018H,3,2021-01-23 15:00:00,34.56478768119395,131.6648721823454,275.6
PL018H,3,2021-01-23 15:17:00,34.57526172144755,131.64681019648799,150.9
PL018H,3,2021-01-23 15:31:00,34.55742359288779,131.64974544980433,28.1
PL017H,3,2021-01-23 16:29:00,34.56290853274202,131.65785829609348,58.1
PL017H,3,2021-01-23 16:47:00,34.56356860883733,131.67800905413438,297.4
PL017H,3,2021-01-23 16:58:00,34.58279819908277,131.66411013071897,175.8
PL018H,4,2021-01-23 17:43:00,34.57104247866242,131.6519106403351,162.2
PL018H,4,2021-01-23 17:56:00,34.54630899468667,131.64520234461372,42.7
PL018H,4,2021-01-23 18:15:00,34.55358860934814,131.6803117712807,286.4
PL017H,4,2021-01-23 19:19:00,34.54611439364768,131.67290361123,337.5
PL017H,4,2021-01-23 19:37:00,34.57080422322992,131.68088004786293,219.2
PL017H,4,2021-01-23 19:48:00,34.55875364470892,131.6523581706941,102.2
PL018H,5,2021-01-23 20:29:00,34.55392045859,131.67484126384787,337.6
PL018H,5,2021-01-23 20:44:00,34.57080171308114,131.67964814102226,221.0
PL018H,5,2021-01-23 20:54:00,34.56382004115043,131.65619269555788,103.2
PL017H,5,2021-01-23 22:01:00,34.57356939688288,131.67052812626355,168.6
PL017H,5,2021-01-23 22:13:00,34.54324094318141,131.66565064177243,30.8
PL017H,5,2021-01-23 22:30:00,34.55594925466217,131.69096983118536,276.7

フィールド調査などでは、測位地点をGPSで落とし、野帳に方位を記入してあとで属性結合してもよいですし、Qfieldなどと連携してデジタル野帳を作っても良いかもしれませんね。

どの処理ツールを使うか？

大変ニッチな分野ではありますが、1990年代後半から様々なVHFテレメトリデータ用のツールが開発されては消えていきました。いくつかのツールについては佐伯・早稲田(2006)にまとまっているので割愛しますが、たとえばAnimal movement(Hooge, P. N. (2000))はArcView3.x上で稼働するエクステンションで、当時わたしと知人でわざわざアラスカ僻地のUSGS研究所まで訪ねていって開発者に「日本語版作らせてくれ」ってお願いしたら「ああいいよ！それで？」と言われたのは良い思い出です()

さて、QGISの拡張機能にはanimal-movementsというタグがあり、

Animove
Radiotrack
がリストアップされてます。

お、これでできるやんといろいろ調べたのですが、Animmoveプラグインに関しては基本的に測位点からその先の解析（縄張りの解析など）を行うためのツールで、今回の地獄の前処理対応にはフィットしません。なお、いろいろ調べると、測位点の深い解析に関してはAnimoveよりもRパッケージのadehabitatHRなどを用いている論文が圧倒的に多い(180/4500ほど)こともあり、adehabitatHRのほうがスタンダードかもしれません、これによるデータ処理にご興味のある方はこちらの先人様の記事などをご参照くださいませ

CF:野生動物の行動圏を推定するパッケージadehabitatHRとGISとの連携

ではRadiotrackプラグインどうかというと、こちらも色々いじってみたものの、もともと発信範囲が短いコウモリなど小型の動物データ処理ということもあり、交差距離が短くないとそもそも三角測量を行わないなどの制限が多いようで、私のデータではどうやってもマチバリ以上の可視化はできませんでした。拡張機能の中身を書き換えるとできそうではありますが、ちょっとめんどくさいかもしれません。

※なお、テレメトリーデータの可視化に特化したTriangulationというプラグインもありましたが、最終更新は2018年9月で、QGIS2.xまでのツールであり書き直すのもちょっとヘビーです。

結局Pythonで処理

手法もデータもニッチであり、なおかつ亜流も多数存在する状況を加味すると、人様のプラグインというパッケージを頼りにするのもなんか違うな、ということでAI氏と~~延々と殴り合い~~入念に議論し結局QGISのPythonコンソールで動くコードを作成しました。
便宜上3本としていますが、このコードを元に改造すればよいかと思います。ただし多角形の重心てどうするんだろう、まあなんとかなるはずです…

また、このコードでもは磁気偏角は考慮していませんでご注意ください。

処理の流れはこんな感じです

コードはこちらになります。初心者向けにコメントは多めにしています。QGIS 3.40.5-Bratislava on MacOSで確認しています。
サクッと見るためのものなのでEPSGは3857に変換してます、そのあたりは適宜変更してください。

# -*- coding: utf-8 -*-
"""
VHF テレメトリ三角測量の可視化スクリプト
機能まとめ：
・観測点（▲）
・方位線
・三角形（交点）
・重心（X）

入力 CSVカラム
name, set_id, datetime, lat, lon, azi_deg
"""

from qgis.core import *
from qgis.PyQt.QtGui import QColor
from PyQt5.QtCore import QVariant
import math, csv, random

# ------------------------------------------------------------
# 入力 CSV のパス（書き換えてください）
# ------------------------------------------------------------
CSV_PATH = "/hoge/hogehoge/hogegegege/vhf_result.csv"

# CSV カラム名
COL_NAME = "name"
COL_ID   = "set_id"
COL_DATETIME = "datetime"
COL_LAT  = "lat"
COL_LON  = "lon"
COL_AZI  = "azi_deg"

REQUIRED_COLS = [COL_NAME, COL_ID, COL_DATETIME, COL_LAT, COL_LON, COL_AZI]

FIXED_LENGTH = 7000  # 方位線の長さ（m）ここも適宜

proj = QgsProject.instance()
crs_wgs  = QgsCoordinateReferenceSystem("EPSG:4326")
crs_proj = QgsCoordinateReferenceSystem("EPSG:3857")
transform_to_proj = QgsCoordinateTransform(crs_wgs, crs_proj, proj)

# ------------------------------------------------------------
# CSV 読み込み（name + set_id で3点をグループ化）
# ------------------------------------------------------------
records = {}

with open(CSV_PATH, newline="", encoding="utf-8") as f:
    reader = csv.DictReader(f)

    for row in reader:
        # 必要なカラムのみ抽出
        filtered_row = {k: row.get(k, "") for k in REQUIRED_COLS}
        
        gid = f"{filtered_row[COL_NAME]}_{filtered_row[COL_ID]}"

        try:
            lat = float(filtered_row[COL_LAT])
            lon = float(filtered_row[COL_LON])
            azi = float(filtered_row[COL_AZI])
        except (ValueError, KeyError) as e:
            print(f"[SKIP] 数値変換エラー: {filtered_row} - {e}")
            continue

        if gid not in records:
            records[gid] = []

        records[gid].append({
            "lat": lat,
            "lon": lon,
            "azi": azi,
            "attributes": filtered_row
        })

print(f"読み込み完了: {len(records)} グループ")

# ------------------------------------------------------------
# レイヤ作成
# ------------------------------------------------------------
obs_layer  = QgsVectorLayer(f"Point?crs={crs_proj.authid()}", "observation_points", "memory")
line_layer = QgsVectorLayer(f"LineString?crs={crs_proj.authid()}", "bearing_lines", "memory")
poly_layer = QgsVectorLayer(f"Polygon?crs={crs_proj.authid()}", "triangles", "memory")
cent_layer = QgsVectorLayer(f"Point?crs={crs_proj.authid()}", "centroids", "memory")

layers = [obs_layer, line_layer, poly_layer, cent_layer]

# 属性追加（必要なフィールドのみ）
fields_to_add = [QgsField(f, QVariant.String) for f in REQUIRED_COLS]
for lyr in layers:
    lyr.dataProvider().addAttributes(fields_to_add)
    lyr.updateFields()

obs_dp  = obs_layer.dataProvider()
line_dp = line_layer.dataProvider()
poly_dp = poly_layer.dataProvider()
cent_dp = cent_layer.dataProvider()

# ------------------------------------------------------------
# 三角測量処理
# ------------------------------------------------------------
triangle_count = 0

for gid, recs in records.items():
    if len(recs) != 3:
        print(f"[SKIP] {gid}: {len(recs)} 点 → 三角形不可")
        continue

    line_geoms = []

    # ------------------------------
    # ① 観測点（▲）＋方位線（→）
    # ------------------------------
    for r in recs:
        pt_wgs = QgsPointXY(r["lon"], r["lat"])
        pt = transform_to_proj.transform(pt_wgs)

        # 観測点
        pf = QgsFeature(obs_layer.fields())
        pf.setGeometry(QgsGeometry.fromPointXY(pt))
        for k, v in r["attributes"].items():
            pf[k] = str(v)
        obs_dp.addFeature(pf)

        # 方位線
        ang = math.radians(r["azi"])
        dx = math.sin(ang) * FIXED_LENGTH
        dy = math.cos(ang) * FIXED_LENGTH
        end_pt = QgsPointXY(pt.x() + dx, pt.y() + dy)

        geom = QgsGeometry.fromPolylineXY([pt, end_pt])
        line_geoms.append(geom)

        lf = QgsFeature(line_layer.fields())
        lf.setGeometry(geom)
        for k, v in r["attributes"].items():
            lf[k] = str(v)
        line_dp.addFeature(lf)

    # ------------------------------
    # ② 方位線の交点 → 三角形
    # ------------------------------
    inter_pts = []
    for g1, g2 in [
        (line_geoms[0], line_geoms[1]),
        (line_geoms[1], line_geoms[2]),
        (line_geoms[2], line_geoms[0]),
    ]:
        inter = g1.intersection(g2)
        if not inter.isEmpty() and inter.type() == QgsWkbTypes.PointGeometry:
            inter_pts.append(inter.asPoint())

    if len(inter_pts) != 3:
        print(f"[WARN] 交点不足: {gid} ({len(inter_pts)}点)")
        continue

    # 三角形作成
    tri = QgsGeometry.fromPolygonXY([inter_pts + [inter_pts[0]]])
    pf = QgsFeature(poly_layer.fields())
    pf.setGeometry(tri)
    for k, v in recs[0]["attributes"].items():
        pf[k] = str(v)
    poly_dp.addFeature(pf)

    # ------------------------------
    # ③ 重心（X）
    # ------------------------------
    cx = sum(p.x() for p in inter_pts) / 3
    cy = sum(p.y() for p in inter_pts) / 3

    cf = QgsFeature(cent_layer.fields())
    cf.setGeometry(QgsGeometry.fromPointXY(QgsPointXY(cx, cy)))
    for k, v in recs[0]["attributes"].items():
        cf[k] = str(v)
    cent_dp.addFeature(cf)

    triangle_count += 1

print(f"三角形生成数: {triangle_count}")

# ------------------------------------------------------------
# 色分け（セットごと＆nameごと）
# ------------------------------------------------------------
def distinct_color(index, total):
    """視認性の高い色を均等に生成（HSV色空間）"""
    hue = (index * 360 / total) % 360
    saturation = 0.7  # 70%
    value = 0.85      # 85%
    
    # HSVをRGBに変換
    import colorsys
    r, g, b = colorsys.hsv_to_rgb(hue / 360, saturation, value)
    return f"#{int(r*255):02x}{int(g*255):02x}{int(b*255):02x}"

# セットごと（name_set_id）に視認性の高いカラー生成（線・測位点・三角形用）
unique_sets = sorted({f"{rec['attributes'][COL_NAME]}_{rec['attributes'][COL_ID]}" 
                      for g in records.values() for rec in g})
set_color_map = {set_id: distinct_color(i, len(unique_sets)) for i, set_id in enumerate(unique_sets)}

# nameごとに視認性の高いカラー生成（重心用）
unique_names = sorted({rec['attributes'][COL_NAME] for g in records.values() for rec in g})
name_color_map = {name: distinct_color(i, len(unique_names)) for i, name in enumerate(unique_names)}

# 観測点 ●（セットごとにカラー）
obs_renderer = QgsCategorizedSymbolRenderer(f"{COL_NAME}||'_'||{COL_ID}", [])
for set_id in unique_sets:
    sym = QgsMarkerSymbol.createSimple({
        "name": "circle", 
        "color": set_color_map[set_id], 
        "size": "3.5",
        "outline_width": "0"
    })
    obs_renderer.addCategory(QgsRendererCategory(set_id, sym, set_id))
obs_layer.setRenderer(obs_renderer)

# 重心 ■（nameごとにカラー・白縁取り）
cent_renderer = QgsCategorizedSymbolRenderer(COL_NAME, [])
for name in unique_names:
    sym = QgsMarkerSymbol.createSimple({
        "name": "square",
        "color": name_color_map[name],
        "size": "5",
        "outline_color": "white",
        "outline_width": "0.8"
    })
    cent_renderer.addCategory(QgsRendererCategory(name, sym, name))
cent_layer.setRenderer(cent_renderer)

# 方位線 →（矢印付き・セットごとにカラー）
line_renderer = QgsCategorizedSymbolRenderer(f"{COL_NAME}||'_'||{COL_ID}", [])
for set_id in unique_sets:
    # ベースラインシンボル
    base = QgsLineSymbol.createSimple({
        "color": set_color_map[set_id], 
        "width": "0.6",
        "capstyle": "flat"
    })
    
    # 矢印マーカー
    arrow = QgsMarkerLineSymbolLayer()
    arrow.setPlacement(QgsMarkerLineSymbolLayer.LastVertex)
    arrow_marker = QgsMarkerSymbol.createSimple({
        "name": "filled_arrowhead",
        "color": set_color_map[set_id],
        "size": "3.5",
        "angle": "0",
        "offset": "0,0"
    })
    arrow.setSubSymbol(arrow_marker)
    
    base.appendSymbolLayer(arrow)
    line_renderer.addCategory(QgsRendererCategory(set_id, base, set_id))
line_layer.setRenderer(line_renderer)

# 三角形 Polygon（セットごとにカラー・薄く）
poly_renderer = QgsCategorizedSymbolRenderer(f"{COL_NAME}||'_'||{COL_ID}", [])
for set_id in unique_sets:
    color = set_color_map[set_id]
    sym = QgsFillSymbol.createSimple({
        "color": color,
        "outline_color": color,
        "outline_width": "0.8",
        "style": "solid"
    })
    # 塗りつぶしの透明度を上げる
    sym.setOpacity(0.3)
    poly_renderer.addCategory(QgsRendererCategory(set_id, sym, set_id))
poly_layer.setRenderer(poly_renderer)

# ------------------------------------------------------------
# set_id のラベル（白縁取り・重心のみ）
# ------------------------------------------------------------
# 重心レイヤーのみラベル表示（[name]-[set]形式）
settings = QgsPalLayerSettings()
settings.fieldName = f'concat("{COL_NAME}", \'-\', "{COL_ID}")'
settings.isExpression = True
settings.enabled = True
settings.placement = QgsPalLayerSettings.AroundPoint

fmt = QgsTextFormat()
fmt.setSize(10)
fmt.setColor(QColor("black"))

buf = QgsTextBufferSettings()
buf.setEnabled(True)
buf.setSize(1.5)
buf.setColor(QColor("white"))
fmt.setBuffer(buf)

settings.setFormat(fmt)

labeling = QgsVectorLayerSimpleLabeling(settings)
cent_layer.setLabeling(labeling)
cent_layer.setLabelsEnabled(True)

# ------------------------------------------------------------
# レイヤ登録
# ------------------------------------------------------------
for lyr in layers:
    lyr.updateExtents()
    proj.addMapLayer(lyr)

print("==============================================")
print("VHF テレメトリ三角測量 完了")
print(f"三角形生成数: {triangle_count}")
print("==============================================")

このコードの最初のほうにある入力CSVのパスをご自分で変更してQGISのPythonコンソールで動かすとこんなデータが作成されます。
個体ID[name]や調査セット[set]ごとに色を変えており、ラベルも表示させています

なお、これらのデータは一時ファイルで出力されますので、QGISを閉じると消えます。処理が終わったら適宜GISファイルにエクスポートをすることをお忘れなく。

今日も頑張る素晴らしきフィールドワーカーたちにも素敵なクリスマスが訪れますように🎅