■概要
「Javaで静岡県ポイントクラウドDBのデータを読み込んで、航空写真と標高PNGを生成する。」で作成した標高PNG画像から、国土地理院で紹介されている「航空レーザ測量データを用いた樹高等のデータ作成」を参考に立木位置や樹冠の検出を試してみます。
■推定DEMの作成
「Javaで静岡県ポイントクラウドDBのデータを読み込んで、航空写真と標高PNGを生成する。」で作成した標高PNGはDSM（デジタル サーフェス モデル）に相当すると思われるので、ここから地形部分のデータを取り出す必要があります。
国土地理院HPの「航空レーザ測量の仕組み」をみると、森林部では、レーザは樹木等で阻害されて地上に届かないこともあり、届くこともあるようです。このため、今回は標高PNGに2ｍメッシュを設定し、各メッシュで最も標高が低い点を地表の標高とみなしました。
以下のようなDemMeshImterpolationProcesserクラスを作成し、DEMに相当する標高PNGを生成しました。
なお、ImageProcessingObservableは進捗状況確認のためのObservableインターフェイスです。
import java.awt.image.BufferedImage;
import java.io.File;
import java.util.Observable;
import javax.imageio.ImageIO;
public class DemMeshImterpolationProcesser extends Observable implements ImageProcessingObservable{
private double[][] dem;
private Cell[][] cell;
private BufferedImage image;
private String message;
public DemMeshImterpolationProcesser(BufferedImage img){
dem=ElevationPngUtil.imgToDouble(img);
image=new BufferedImage(img.getWidth(),img.getHeight(),BufferedImage.TYPE_INT_RGB);
}
public BufferedImage getImage(){
return image;
}
public void drawMinH(){
for(int i=0;i<cell.length;i++){
for(int j=0;j<cell[i].length;j++){
cell[i][j].setRGB(ElevationPngUtil.getRGB(cell[i][j].getMinH()));
}
}
}
public void createImageCell(int size){
int ww=dem.length;
int hh=dem[0].length;
cell=new Cell[ww/size][hh/size];
int x=0;
int y=0;
for(int i=0;i<cell.length;i++){
for(int j=0;j<cell[i].length;j++){
cell[i][j]=new Cell();
cell[i][j].x=new int[size];
cell[i][j].y=new int[size];
for(int m=0;m<size;m++){
cell[i][j].y[m]=y;
y=(y+1)%hh;
}
for(int m=0;m<size;m++){
cell[i][j].x[m]=x+m;
}
if(y==0)x=x+size;
}
}
}
public void output(File f)throws Exception{
ImageIO.write(image, "png", f);
update("出力終了");
}
protected void update(String mes){
message=mes;
setChanged();
notifyObservers();
}
@Override
public String progress() {
return message;
}
protected class Cell{
int[] x;
int[] y;
double getMinH(){
double min=Double.MAX_VALUE;
for(int i=0;i<x.length;i++){
for(int j=0;j<y.length;j++){
if(Double.isNaN(dem[x[i]][y[j]])||dem[x[i]][y[j]]<0)continue;
min=Math.min(min, dem[x[i]][y[j]]);
}
}
if(min==Double.MAX_VALUE){
return 0;
}else{
return min;
}
}
void setRGB(int rgb){
for(int i=0;i<x.length;i++){
for(int j=0;j<y.length;j++){
image.setRGB(x[i], y[j], rgb);
}
}
}
}
}
上記で生成したDEM標高PNGです。（以下はサイズを縮小したJPGです。）
なお、標高PNG生成後、5×5のガルシアンフィルタで平滑化しました。
■推定DSMの作成
前回作成した標高PNGはデータの無いピクセルも多くあるので、TINを生成して補間しました。
TINの作成及び補間には「Tinfour」を利用しました。
「Tinfour」のライセンスは「the Apache License, Version 2.0」です。
/**
* TIN生成
*/
public void process(){
for(int i=0;i<img.getWidth();i++){
for(int j=0;j<img.getHeight();j++){
int rgb=img.getRGB(i, j);
double hh=ElevationPngUtil.getZ(rgb);
if(Double.isNaN(hh))continue;
Vertex v=new Vertex(i,j,hh,index++);
list.add(v);
}
}
tin=new IncrementalTin();
tin.add(list, null);
}
/**
* TINによる補間
* @return
*/
public BufferedImage interpolation(){
BufferedImage ret=new BufferedImage(img.getWidth(),img.getHeight(),BufferedImage.TYPE_INT_RGB);
setInitVal(ret);
TriangularFacetInterpolator tfi=new TriangularFacetInterpolator(tin);
VertexValuatorDefault vvd=new VertexValuatorDefault();
double n=img.getWidth()*img.getHeight();
for(int i=0;i<img.getWidth();i++){
for(int j=0;j<img.getHeight();j++){
double hh=tfi.interpolate(i, j, vvd);
if(hh>0)ret.setRGB(i, j, ElevationPngUtil.getRGB(hh));
}
}
return ret;
}
TINで内層したDSM標高PNGは以下のとおりです。（以下はサイズを縮小したJPGです。）
■推定DCHMの作成
DSMからDEMを差し引いたのが、樹木等の高さの標高モデル（DCHM）となります。
DSM標高データからDEM標高データを差し引いて作成したDCHM標高PNGは以下のとおりです。（以下はサイズを縮小したJPGです。）
なお、DCHM標高PNGから作成した標高図（樹高図）は以下のとおりです。
■樹木ピークの検出
「航空レーザ測量の仕組み」では3ｍメッシュで最大値の地点を検出する方法が紹介されています。
今回は半径2.0ｍの探索円を設定し、以下のコードで頂点位置を検出しました。
private void findPeaks(){
for(double z=highH;z>=lowH;z=z-1){
double dv=(++vv)/nn*100;
for(int i=0;i<width;i++){
for(int j=0;j<height;j++){
if(ck[i][j]!=0)continue;
if(val[i][j]>=z){
Ellipse2D e=new Ellipse2D.Double(i-this.radius, j-this.radius, this.radius*2, this.radius*2);
if(isHiMax(e,val[i][j])){
setPeak(e);
points.add(new Point3d(e.getCenterX(),e.getCenterY(),val[i][j]));
}
}
}
}
}
}
private boolean isHiMax(Ellipse2D e,double v){
int xs=(int)Math.max(e.getX(), 0);
int ys=(int)Math.max(e.getY(), 0);
int xe=(int )Math.min(e.getX()+e.getWidth(), width);
int ye=(int)Math.min(e.getY()+e.getHeight(),height);
for(int i=xs;i<xe;i++){
for(int j=ys;j<ye;j++){
if(val[i][j]>v)return false;
}
}
return true;
}
private void setPeak(Ellipse2D e){
int xs=(int)Math.max(e.getX(), 0);
int ys=(int)Math.max(e.getY(), 0);
int xe=(int )Math.min(e.getX()+e.getWidth(), width);
int ye=(int)Math.min(e.getY()+e.getHeight(),height);
for(int i=xs;i<xe;i++){
for(int j=ys;j<ye;j++){
if(ck[i][j]!=0)continue;
if(e.contains(i, j))ck[i][j]=1;
}
}
}
上記で検出したHCHMピーク点（立木位置？）は以下のとおりです。
画像の範囲（674m×600ｍ）で12392点のピーク点（立木位置？）が検出されました。
■樹冠の推定
HCHMピーク点を樹木の位置として、「Tinfour」を用いてボロノイ分割を実施しました。
public List<ThiessenPolygon> process(){
BoundedVoronoiBuildOptions options = new BoundedVoronoiBuildOptions();
options.enableAutomaticColorAssignment(true);
BoundedVoronoiDiagram diagram = new BoundedVoronoiDiagram(vertList, options);
List<ThiessenPolygon> polyList=diagram.getPolygons();
return polyList;
}
以下のボロノイ分割結果を示します。
拡大すると、写真で確認できる樹冠の影とピーク・領域がだいたい良い感じ分割されているように見えます。
■ピーク点の確認
ちゃんとピークが検出されているか等を確認するため、樹冠と推定したDCHMピクセル情報をJyz3dで可視化してみました。
何か所か確認しましたが、だいたい下図のような形でピークが検出されていました。
■最後に
静岡県ポイントクラウドDBデータからJavaで樹高等の解析を試してみましたが、なんとなく、それっぽく検出できました。
検証データがないので、遊びの範囲ではありますが、広葉樹等では厳しいですが、スギ植林などではそれなりに樹木本数や樹冠を検出できるように感じました。
なお、「航空レーザ測量データを用いた樹高等のデータ作成」では『航空レーザ測量による樹高データは、点群密度が高い常緑針葉樹林のデータを除き、一般的には、個々の樹木の実測値とはあまり合いません。』と述べられており、また、『しかし、30m四方程度の範囲で平均値を取ると、相関が良くなります。』と述べられています。
静岡県ポイントクラウドDBのように、測量や調査の結果がもっと広く公開されれば、いろんな事に利用できそうです。