Three.jsで樹頂点ポイントデータをもとに森林を3Dビジュアル化してみる

demoページ

はじめに

鳥取県が森林情報に関するデータを公開していました。

（本データは「Tottori Forestory Innovation+（TFI+）」において参加者へ提供するデータのサンプルデータです。）

とのことで、どうやら森林・林業分野においての開発コンペ用のサンプルデータらしいですが、このデータの中に樹頂点ポイントデータがありました。

QGISで表示した樹頂点ポイントデータ（ヒノキ、スギ、マツ）

鳥取県の一部の山林のみのデータになりますが、こういった樹木一本一本がポイントデータとなって公開されているのはなかなかないので、せっかくなのでこのデータを使って遊んでみました。

樹頂点ポイントデータは、森林内の個々の樹木の最も高い位置、すなわち樹頂を示す地点のデータです。LiDARによるレーザ計測で得た樹冠高データ（DCHM）から作られてるポイントデータなので、樹木に関する６つの項目（樹冠投影面積、樹冠表面積、樹冠体積、樹高、樹冠長、樹冠長率）を含んでいます。もちろん位置情報もあります。

今回は、このポイントデータからThree.jsで3Dの立木にしちゃいます。といっても単純な3D図形による描画になります。

データの加工

今回使用するデータはそのままではちょっと扱いにくいので、コマンドラインでデータ加工していきます。
（本題のthree.jsの部分は後半になるので、データ加工に興味ない方は読み飛ばしてください。）

今回使用するのは以下のファイルです。

• ヒノキ_樹頂点_05LE904.shp
• スギ_樹頂点_05LE904.shp
• マツ_樹頂点_05LE904.shp
• 05LE904.tif（DEM）
• 05LE904.tif（微地形表現図）

（同じファイル名がありますがそれぞれ違うものです。ややこしいので加工時にファイル名も分けます。）

データには加工にはGDALとcsvtojsonを使用してコマンドラインで加工しますが、csvtojsonは最後の圧縮したjsonに出力した時のみ使用したのでほぼほぼGDALによる加工になります。
また、なるべくJavaScriptのコードの中に計算式を書くのをなるべく省きたかったので、データ加工の段階で必要な情報を計算してデータに書き込むことにしました。

樹頂点ポイントデータの加工（shp→JSON）

データの属性テーブル（スギ）

CRSはEPSG:2447

まずはこの３樹種の樹頂点データを一つにまとめちゃいます。

SHPファイルのデータをWEBに扱いやすいデータに加工します。単純にGeoJsonにしてしまえばJavaScript側で普通に読めますが、今回は地理情報データである必要はないので、コンパクトなJSONにしちゃいます。

まずは、それぞれの樹頂点データを一つのshpファイルにまとめます。ogrmerge.pyコマンドでヒノキ_樹頂点_05LE904.shp、マツ_樹頂点_05LE904.shp、スギ_樹頂点_05LE904.shpをマージして、新たにshpファイルを作成します。この時に-lco ENCODING=Shift_JISでShift_JISエンコーディングを指定する必要がありました。

ogrmerge.py -single -o treetop.shp -nln treetop -lco ENCODING=Shift_JIS ヒノキ_樹頂点_05LE904.shp マツ_樹頂点_05LE904.shp スギ_樹頂点_05LE904.shp

次に今回の3D可視化に必要な情報を属性テーブルに書き込みます。まずは、胸高直径の値を求めます。

胸高直径（きょうこうちょっけい、dbh: diameter at breast height）とは、樹木の幹の太さを測る指標の一つで、樹幹の直径を胸高1.2mの高さで測定した値のことを指します。樹木の大きさを比較しやすくするために用いられます。

胸高直径は（ 形状比 / 100）/ 樹高 で計算できます。

参考資料

形状比とは樹高（m）を直径（cm）で割り100を掛けたものです。その値が大きいほどヒョロヒョロな状態であることを表します。形状比が高い木は、風雪害に弱いと言われています（100を超えるような木はかなり危険）。

もともと胸高直径はセンチ単位で表すのが主流ですが、他の単位と合わせるためにm単位の数値に修正します。

ogrinfoコマンドでSQL文を使用して胸高直径のフィールドをALTER TABLE文でADD COLUMNにより追加します。

ogrinfo -sql "ALTER TABLE treetop ADD COLUMN 胸高直径 VARCHAR" treetop.shp

そしたら-dialect sqliteオプションを付けて、UPDATE文で形状比のテーブルの値から計算して胸高直径の値をSETで追加します。ROUND関数により小数点を第3まで丸めています。

ogrinfo -dialect sqlite -sql "UPDATE treetop SET 胸高直径 = ROUND(樹高 / ( 形状比 / 100 ) / 100, 3)" treetop.shp

次に、樹冠半径の情報を書き込むためのフィールドを追加します。

「樹冠」とは、樹木の上部にある枝と葉の集まりを指します。「樹冠半径」は、この樹冠の中心からその外縁までの距離を指します。平面的な視点から見たとき、樹冠の形が円や楕円に近い場合、その半径を「樹冠半径」として考えることができます。

樹冠半径の値を書き込みます。樹冠半径は√樹冠投影面積 / πで求められます。
ポイントデータに含まれる「樹冠面積」のフィールドが樹冠投影面積の値だったのでここから計算します。
（別で「樹冠表面積」という似たような名前のフィールドが含まれてましたが、樹冠投影面積を「樹冠体積」のフィールドの値から算出したところ、「樹冠面積」の値と大体一致していたので、こちらが樹冠投影面積であっていると思います。）

樹冠半径のフィールドを追加し、SQRT( 樹冠面積 / PI() )で計算してROUND関数で小数点を丸めて計算結果の値を樹冠半径のテーブルに追加します。

ogrinfo -sql "ALTER TABLE treetop ADD COLUMN 樹冠半径 VARCHAR" treetop.shp

ogrinfo -dialect sqlite -sql "UPDATE treetop SET 樹冠半径 = ROUND( SQRT( 樹冠面積 / PI() ), 1)" treetop.shp

次にogr2ogrコマンドでshpをgeojsonに変換します。元でーたの座標値に余分な0が含んでおり、無駄に冗長な精度値になっていたので-lco SIGNIFICANT_FIGURES=8で、緯度経度の値を丸めてることで、余分な数値を削除します。また、マルチポイントのジオメトリになっていて使いづらかったので、-nlt POINT オプションでポイントのジオメトリに変換します。さらにSQL文のSELECT構文で必要な属性フィールドだけを抜き出して、フィールド名を短いアルファベットに省略することで容量を軽くします。

ogr2ogr -f GeoJSON -nlt POINT -lco SIGNIFICANT_FIGURES=8 treetop.geojson treetop.shp -dialect sqlite -sql "SELECT 中樹種ID AS ID, 樹高 AS H, 樹冠長 AS Cl, 胸高直径 AS D, 樹冠半径 AS Bl, geometry FROM treetop"

抜き出したフィールド名と、書き換えた省略アルファベット

Original	Alias
中樹種ID	ID
樹高	H
樹冠長	Cl
胸高直径	D
樹冠半径	Bl

次に、このgeojsonをCSVに変換し、GEOMETRY=AS_XYのオプションによりポイントの位置情報をテーブルに追加します。

ogr2ogr -f CSV treetop.csv treetop.geojson -lco GEOMETRY=AS_XY

ここでGDALの出番は終わりで最後にcsvtojsonで処理します

CSVをcsvtojsonコマンドでJSONに変換します。JSONのvalueの型が文字列型にならないように--checkType=trueをつけてます。

csvtojson treetop.csv --checkType=true > treetop.json

これで一つのJSONにまとめた樹頂点ポイントデータができました。もともと3つ合わせて15MBくらいでしたが3.1MBまで圧縮できました。。

JSONの中身はこんな感じになってます。位置情報と先ほど抜き出したフィールドの情報が入ってます。

[
{"X":2053.75,"Y":-59249.25,"ID":2,"H":13.7,"Cl":7.9,"D":0.225,"Bl":2},
{"X":2459.75,"Y":-59249.25,"ID":2,"H":10.7,"Cl":2.8,"D":0.168,"Bl":1.7},
{"X":2469.75,"Y":-59249.25,"ID":2,"H":10.7,"Cl":5.5,"D":0.157,"Bl":1.5},
{"X":2475.25,"Y":-59249.25,"ID":2,"H":8.8,"Cl":3.7,"D":0.147,"Bl":1.6},
{"X":2691.75,"Y":-59249.25,"ID":2,"H":22.8,"Cl":13.5,"D":0.43,"Bl":2.9},
{"X":2645.25,"Y":-59248.75,"ID":2,"H":21.3,"Cl":10.4,"D":0.359,"Bl":2.5},
{"X":2480.75,"Y":-59247.75,"ID":2,"H":9.6,"Cl":5.4,"D":0.149,"Bl":1.6},
{"X":2671.75,"Y":-59247.75,"ID":2,"H":18.6,"Cl":6.7,"D":0.316,"Bl":2.4},
{"X":2447.75,"Y":-59247.25,"ID":2,"H":9.3,"Cl":5.4,"D":0.208,"Bl":2.2},
{"X":2468.75,"Y":-59247.25,"ID":2,"H":10.6,"Cl":6.1,"D":0.178,"Bl":1.8},

/// 以下省略 ///

shp→geojson→csv→jsonというわけわからん変換手順を踏んでいますが、GDALのそれぞれのドライバー固有のオプションを使用してジオメトリの情報とテーブルの情報をうまく結合したかったのでこのような手順になりました。

DEMデータの加工

QGISで表示したDEM

詳細情報

なぜか座標系（CRS）の情報が入ってませんが、樹頂点のポイントの座標系と同じEPSG:2447で間違いないでしょう。0.5mメッシュで解像度のいいデータなのですが、このまま使用するとブラウザがクラッシュしたため、今回はgdalwarpコマンドで解像度を5mに落とします。

gdalwarp -tr 5 5 05LE904.tif dem_5m.tif

微地形表現図データの加工(Geotiff→webp)

QGISで表示した微地形表現図

微地形表現図のgeotiffを加工します。これはthree.jsでオブジェクトのテクスチャ用に使用するのでgdal_translateコマンドを使ってwebpに変換します。

gdal_translate 05LE904.tif microtopographic.webp

データ加工はこれで終わりです。

Three.jsに表示する

ここからはThree.jsで先ほどのデータを使って実装していきます。まずはレンダリングに必要なキャンバスを用意し、シーンを作成して、オブシェクトの表示に必要なカメラ、ライト、コントロールを追加します。

import * as THREE from 'three';
import { OrbitControls } from 'three/examples/jsm/controls/OrbitControls.js'; 

// 画面サイズ
const sizes = {
    width: window.innerWidth,
    height: window.innerHeight,
};

// キャンバスを取得
const canvas = document.querySelector('.webgl');

// シーンの作成
const scene = new THREE.Scene();

// ライト（３つ）
const TopDirectionalLight = new THREE.DirectionalLight(0xffffff, 0.5);
TopDirectionalLight.position.set(0, 1000, 0);

const SideDirectionalLight = new THREE.DirectionalLight(0xffffff, 0.5);
SideDirectionalLight.position.set(0, 0, 1000);

const ambientLight = new THREE.AmbientLight(0xffffff, 0.5);
scene.add(TopDirectionalLight, SideDirectionalLight, ambientLight);

// カメラ
const camera = new THREE.PerspectiveCamera(75, sizes.width / sizes.height, 0.1, 10000);

// 極半径 r
const r = 1500;

// 方位角
const azimuthAngle = THREE.MathUtils.degToRad(-62);

// 仰角
const elevationAngle = THREE.MathUtils.degToRad(65);

// x, y, z座標を計算
const x = r * Math.sin(elevationAngle) * Math.cos(azimuthAngle);
const y = r * Math.cos(elevationAngle);
const z = r * Math.sin(elevationAngle) * Math.sin(azimuthAngle);

// カメラの位置を設定
camera.position.set(x, y, z);

// コントロール
const controls = new OrbitControls(camera, canvas);
controls.enableDamping = true;

// レンダー
const renderer = new THREE.WebGLRenderer({
    canvas: canvas,
});
renderer.setSize(sizes.width, sizes.height);
renderer.setPixelRatio(Math.min(window.devicePixelRatio, 2));

// アニメーション
const animate = () => {
    requestAnimationFrame(animate);

    controls.update();
    renderer.render(scene, camera);
};
animate();

// 画面リサイズ時にキャンバスもリサイズ
window.addEventListener(
    'resize',
    () => {
        camera.aspect = window.innerWidth / window.innerHeight;
        camera.updateProjectionMatrix();
        renderer.setSize(window.innerWidth, window.innerHeight);
    },
    false,
);

背景は、いい感じにしたいので、Skyオブジェクトでシーンの中に空を作成します。
この時に、太陽の方向に合わせて、先ほどシーンに追加したライト一つを移動させてます。

import { Sky } from 'three/examples/jsm/objects/Sky.js';

// Skyオブジェクトを作成
const sky = new Sky();

// Skyのスケールを設定
sky.scale.setScalar(450000);

// Skyオブジェクトをシーンに追加
scene.add(sky);

// Skyオブジェクトのマテリアルのuniformsを取得
const uniforms = sky.material.uniforms;

// 大気の透明度
uniforms['turbidity'].value = 10;

// 空の青さの度合い
uniforms['rayleigh'].value = 3;

// 太陽光の散乱度
uniforms['mieCoefficient'].value = 0.005;
uniforms['mieDirectionalG'].value = 0.7;

// 太陽の位置や角度を制御するためのパラメータ
const parameters = {
    inclination: 0.49, // 太陽の傾斜角
    azimuth: -32.4, // 太陽の方位角
    elevation: 2, // 太陽の高度（地平線からの角度）
};

// 太陽の位置を表すベクトルの初期化
const sun = new THREE.Vector3();

// 球座標系を使用して太陽の位置を決定
const phi = THREE.MathUtils.degToRad(90 - parameters.elevation);
const theta = THREE.MathUtils.degToRad(parameters.azimuth);
sun.setFromSphericalCoords(1, phi, theta);

// 計算された太陽の位置をuniformに設定
uniforms['sunPosition'].value.copy(sun);

// DirectionalLightの位置を太陽の位置と一致するように設定
SideDirectionalLight.position.copy(sun);

できたシーン

このシーンに地形をつくります。

DEMから地形オブシェクトを作成する

DEMを読み込むためにgeotiff.jsを使用します

npm install geotiff

import { fromArrayBuffer } from 'geotiff';

// ラスター処理
const loadRasterData = async () => {
    // ラスターデータの読み込み
    const filepath = '../dem_5m.tif';
    const response = await fetch(filepath);
    const arrayBuffer = await response.arrayBuffer();
    const tiff = await fromArrayBuffer(arrayBuffer);
    const image = await tiff.getImage();

    // ラスターデータを取得
    const rasters = await image.readRasters();

    // ピクセル解像度を取得
    const modelPixelScale = await image.getFileDirectory().ModelPixelScale;
    const dx = modelPixelScale[0];
    const dy = modelPixelScale[1];

    // ラスターの高さと幅を取得
    const tiffWidth = image.getWidth();
    const tiffHeight = image.getHeight();

    // ピクセル幅の取得
    const pixelWidth = image.getFileDirectory().ModelPixelScale[0];
    const pixelHeight = image.getFileDirectory().ModelPixelScale[1];

    // 微地形表現図のテクスチャ
    const texture = new THREE.TextureLoader().load('../microtopographic.webp');

    // マテリアルの作成
    const planematerial = new THREE.MeshLambertMaterial({
        side: THREE.DoubleSide, // 裏面もテクスチャを表示
        map: texture,
    });

    // PlaneGeometryの作成。ラスターのピクセル数分の頂点を持つ
    const planegeometry = new THREE.PlaneGeometry(tiffWidth * pixelWidth, tiffHeight * pixelHeight, tiffWidth - 1, tiffHeight - 1);

    // デフォルトは縦向きなのでをX軸を中心に90度に回転
    planegeometry.rotateX(-Math.PI / 2);

    // メッシュの作成
    const planegmesh = new THREE.Mesh(planegeometry, planematerial);

    //PlaneGeometryの頂点座標をDEMの値に置き換える
    const vertices = planegeometry.attributes.position.array;
    const demData = rasters[0];

    // DEMの値の最小値を計算
    const minValue = demData.reduce((min, value) => Math.min(min, value), Infinity);

    // DEMの値を元に、vertices配列のy座標を更新して頂点を立ち上げる
    let j = 0;
    demData.forEach((value, i) => {
        vertices[j + 1] = value;
        j += 3;
    });

    // 立ち上げたPlaneGeometryの底面が原点0になるようにジオメトリを下げる
    planegmesh.position.y = -minValue;

    // シーンにPlaneGeometryを追加
    planegeometry.name = 'dem';
    scene.add(planegmesh);
};

// ラスターの読み込み処理を実行
loadRasterData().catch((error) => {
    console.error('Error loading raster data:', error);
});

解説

geotiff.jsでDEMのラスターを読み込み、ラスターのサイズ（幅、高さ）とピクセルサイズ、ピクセル数を取得します。
地形を作るための土台になるPlaneGeometryを作成します。このジオメトリのサイズはラスターで取得したサイズを割り当て、ジオメトリの中の頂点数はDEMのピクセル数から１を差し引いた数にしています。これは、内側の頂点数よりもピクセル数の方が1つ分多くなるためです。

そして、このジオメトリの各頂点のZ軸（高さ）にDEMの値を割り当てて頂点を立ち上げることで地形の起伏を再現します。こうすることで、Three.js上でDEMから地形オブジェクトを作成することができます。

そして、地形オブジェクトのテクスチャには、先ほど加工した微地形表現図の画像に設定してます。

ということで地形オブジェクトができました。

参考

樹頂点データから樹木オブジェクトの作成

次に、この地形に樹木を生やします。樹頂点データのポイントは約40000ポイントあるのでおよそ40000本の樹木オブシェクトを作成する必要があります。

大量のオブジェクトを作成するため、オブジェクトのマージ処理用にこちらのBufferGeometryUtilsアドオンをディレクトリに追加しておきます。

樹木オブジェクトを地形オブジェクトの上に生やすために、レイキャスティング処理を大量に行うため、レイキャスティングを高速化するモジュールをインストールします

npm i three-mesh-bvh

モジュールをインポートをし、BVH(Bounding Volume Hierarchy)を有効にしておきます。

import * as BufferGeometryUtils from 'three/examples/jsm/utils/BufferGeometryUtils.js'
import { MeshBVH, acceleratedRaycast } from 'three-mesh-bvh';

// BVHの高速化されたレイキャストを有効にする
THREE.Mesh.prototype.raycast = acceleratedRaycast;

先ほどのラスター処理の関数の中に以下を記述します。

// ラスター処理
const loadRasterData = async () => {

    /// 省略 ///

    // planegeometryに対してBVH (Bounding Volume Hierarchy) を作成
    planegeometry.boundsTree = new MeshBVH(planegeometry);

    // ラスターの中心座標を取得（別の処理で使用）
    const [ulx, uly] = image.getOrigin();
    const centerX = ulx + (tiffWidth / 2) * dx;
    const centerY = uly - (tiffHeight / 2) * dy;
    const center = [centerX, centerY];

    // DEMの値の最大値を計算
    const maxValue = demData.reduce((max, value) => Math.max(max, value), -Infinity);

    // レイキャスティングの原点のY軸を計算（別の処理で使用）
    const rayOriginY = maxValue - minValue;

    // 樹木オブジェクトを作成
    createTreeGeometry(planegmesh, center, rayOriginY);

};

地形ジオメトリ（planegeometry）に対してBVHを作成し、boundsTreeプロパティにセットすることで、ジオメトリの衝突検出やレイキャスティングのパフォーマンスを向上させています。

また、後述するcreateTreeGeometry関数に３つの引数を渡してます

planegmesh
地形オブジェクトに対してレイキャスティングするために使用します。

center
ラスターの中心の位置座標です。樹頂点ポイントをワールド座標に合わせるためのオフセットに使用します。

rayOriginY
レイキャスティングの起点の高さ（Y）です。地形オブシェクトの各頂点から最大の高さの値を渡してます。

createTreeGeometry関数の処理を書いていき、樹木オブジェクトを作成してシーンに追加します。

const createTreeGeometry = (planegmesh, center, rayOriginY) => {

    // Raycasterの設定
    const raycaster = new THREE.Raycaster();
    raycaster.firstHitOnly = true;

    // ファイルloaderをインスタンス化
    const loader = new THREE.FileLoader();

    // 樹木のデータを格納する配列
    const treetopData = [];

    loader.load(
        '../treetop.json',
        (data) => {
            JSON.parse(data).forEach((point) => {
                // 座標
                const getPoint = [point.X - center[0], point.Y - center[1]];
                // 樹種(中樹種ID)
                const treeSpecies = point.ID;
                // 枝下高
                const trunkHeight = point.H - point.Cl;
                // 幹半径
                const trunkRadius = point.D / 2;
                // 樹冠長
                const crownLength = point.Cl;
                // 樹冠半径
                const crownRadius = point.Bl;

                // 幹のジオメトリ
                const trunkGeometry = new THREE.CylinderGeometry(trunkRadius, trunkRadius, trunkHeight, 6);
                // 樹冠のジオメトリ
                const crownGeometry = new THREE.ConeGeometry(crownRadius, crownLength, 6);

                // レイキャスターの原点と方向を設定
                const rayDirection = new THREE.Vector3(getPoint[0], -1, -getPoint[1]);
                const rayOrigin = new THREE.Vector3(getPoint[0], rayOriginY, -getPoint[1]);
                raycaster.set(rayOrigin, rayDirection.sub(rayOrigin).normalize());

                //中樹種IDによって樹冠の色を変更
                let crownColor;
                switch (treeSpecies) {
                    case 1: // スギ
                        crownColor = 'rgb(74, 230, 77)';
                        break;
                    case 2: // ヒノキ
                        crownColor = 'rgb(0, 122, 6)';
                        break;
                    case 3: // マツ
                        crownColor = 'rgb(33, 164, 114)';
                        break;
                    default:
                        break;
                }
                
                
                const geometries = [
                    { geometry: trunkGeometry, defaultY: 0, height: trunkHeight / 2, color: 'rgb(232, 132, 1)' },
                    { geometry: crownGeometry, defaultY: 0, height: crownLength / 2 + trunkHeight, color: crownColor },
                ];

                geometries.forEach((geoData) => {

                    // PlaneGeometry (planegmesh) との交差を検出
                    const intersects = raycaster.intersectObject(planegmesh);
                    
                    // 交差が検出された場合、ジオメトリを交差点のy座標に合わせる
                    if (intersects.length > 0) {
                        geoData.geometry.translate(getPoint[0], intersects[0].point.y + geoData.height, -getPoint[1]);
                    } else {
                        geoData.geometry.translate(getPoint[0], geoData.defaultY, -getPoint[1]);
                    }
                    
                    // 与えられた色情報を元に、各頂点に色を設定
                    const color = new THREE.Color(geoData.color);
                    const colorsArray = Array.from({ length: geoData.geometry.getAttribute('position').count }).flatMap(() => [color.r, color.g, color.b]);
                    geoData.geometry.setAttribute('color', new THREE.Float32BufferAttribute(colorsArray, 3));

                    treetopData.push(geoData.geometry);
                });
            });

            // 樹木のマテリアル
            const treeMaterial = new THREE.MeshLambertMaterial({
                vertexColors: true,
            });

            // ジオメトリをマージ
            const treeSphere = new THREE.Mesh(BufferGeometryUtils.mergeGeometries(treetopData), treeMaterial);

            treeSphere.name = 'treeSphere';
            scene.add(treeSphere);
        },

        (err) => {
            console.error(err);
        },
    );
};

解説

加工した樹頂点のJSONファイルをFileLoaderで読み込みます。JSONには、樹木の各特性（位置、種類、高さ、直径など）に関する情報が含まれています。その情報に基づいて樹木の3Dジオメトリ（形状）を生成しています。

樹木を表現するジオメトリは、幹部分をCylinderGeometry。樹幹部分をConeGeometryで作成します。頂点数はなるべく少量にしています。また、樹種によってConeGeometryに別々の頂点カラーを付与して、樹冠を色分けしてます。

樹木オブジェクトの構造

樹木オブジェクトに割り当てる値

樹木の位置情報は平面直角座標系（EPSG:2447）の座標値を含んでいるため、ワールド座標の原点0から離れた位置になります。ですが、すでにワールド座標系上に平面直角座標系のスケールでプロットされているDEM（地形オブジェクト）の中心点はワールド座標系の0の位置になっています。このDEMの元々の中心の座標値（先ほど渡した変数center）を平面直角座標系の樹木のポイントの座標値から減算することにより、ワールド座標上にプロットされているDEMとの座標値がぴったり重なるようになります。
そして、樹木の根本の高さの位置（Y座標）はレイキャスティングにより、planegmesh（地形オブジェクト）の衝突検知によって求めます。

樹木オブジェクトのY座標値について

最後に、全ての樹木ジオメトリをマージして、一つのジオメトリとして扱うことで描画のパフォーマンスを向上させてます。

シーンに追加した樹木オブジェクト

実際の計測データをもとに生成したので、木一本一本の形状が異なります（データがおかしいものが混じっていて異形な樹木もあったりします）。

樹木の下にカメラを潜り込ませた様子。
幹の太さの違いはちょっと分かりづらいので、数値を乗算して差を大きくすれば分かりやすくなると思います。

クリックして伐採する機能とかつけて、間伐シュミレーションとかできそう・・・。

おわりに

今回はジオメトリの大量生成やオブジェクトの頂点座標の操作をCPU処理で実装してるので、最初の描画に少し時間がかかります。この辺りの部分はカスタムシェーダーを利用してGPUでの処理に変更すれば、表示の高速化が期待できるかもしれません。