Javaで静岡県ポイントクラウドDBのデータを読み込んで、航空写真と標高PNGを生成する。

## ■概要
最近、仕事半分、趣味半分で地理情報処理の勉強をしています。
その中でオープンな点群データがないか検索していたら、[静岡県ポイントクラウドDB](https://pointcloud.pref.shizuoka.jp/)に行きつきました。
[静岡県ポイントクラウドDB](https://pointcloud.pref.shizuoka.jp/)では、Lasデータセットでデータが公開されているので、JavaでLasを読み込んで航空写真と標高データ画像（標高PNG）を生成してみました。

## ■標高PNG処理
標高のデータ化には[国土地理院標高タイルで使用されている標高PNG](https://maps.gsi.go.jp/development/demtile.html)を用いました。
標高PNGの論文はこちらの様です。
→「[PNG標高タイル －Web利用に適した標高ファイルフォーマットの考察と実装－（西岡芳晴・長津樹理、情報地質第26巻第4号、2015）](https://www.jstage.jst.go.jp/article/geoinformatics/26/4/26_155/_pdf)」

地理院HP及び論文を参考に、ElevationPngUtilクラスを作成しました。

~~~java:ElevationPngUtil

import java.awt.Color;
import java.awt.geom.AffineTransform;
import java.awt.geom.Point2D;
import java.awt.image.BufferedImage;
import java.io.BufferedWriter;
import java.io.File;
import java.io.FileWriter;
import java.io.IOException;

/**
 * 標高PNGユーティリティクラス
 */
public class ElevationPngUtil {

	private static final int P8=256;
	private static final int P16=65536;
	private static final int P23=8388608;
	private static final int P24=16777216;
	private static final double U=0.01;

	private ElevationPngUtil(){}

	/* 標高PNG NA値*/
	public static int NA=P23;

	/* RGB値をIntに変換*/
	private static int rgb2Int(int[] c){
		return new Color(c[0],c[1],c[2]).getRGB();
	}

	/**
	 * 標高値を標高PNG画素値に変換
	 * @param z 標高値（ｍ）
	 * @return 標高PNG画素値
	 */
	public static int getRGB(double z){
		if(Double.isNaN(z))return P23;
		return rgb2Int(getRGBColor(z));
	}

	/* 標高値をRGB:int[]に変換 */
	private static int[] getRGBColor(double z){
		if(z<=0)return new int[]{128,0,0};
		int i=(int)Math.round(z*100);
		int r=i >> 16;
		int g=i-(r << 16) >> 8;
		int b=i-((r << 16)+(g << 8));
		return new int[]{r,g,b};
	}

	/**
	 * 標高PNG画素値を標高値に変換
	 * @param intColor 標高PNG画素値
	 * @return 標高値（ｍ）
	 */
	public static double getZ(int intColor){
		Color c=new Color(intColor);
		int r=c.getRed();
		int g=c.getGreen();
		int b=c.getBlue();
		int x=r*P16+g*P8+b;
		if(x<P23){
			return U*(double)x;
		}else if(x>P23){
			return U*(double)(x-P24);
		}else{
			return Double.NaN;
		}
	}

	public static double[] getMinMaxZ(BufferedImage png){
		double min=Double.MAX_VALUE;
		double max=-Double.MAX_VALUE;
		for(int i=0;i<png.getWidth();i++){
			for(int j=0;j<png.getHeight();j++){
				double h=getZ(png.getRGB(i, j));
				if(Double.isNaN(h))continue;
				min=Math.min(h, min);
				max=Math.max(h, max);

			}
		}
		return new double[]{min,max};
	}

	/**
	 * 標高PNGを標高の配列（double[][]）に変換
	 * @param png
	 * @return
	 */
	public static double[][] imgToDouble(BufferedImage png){
		double[][] ret=new double[png.getWidth()][png.getHeight()];
		for(int i=0;i<ret.length;i++){
			for(int j=0;j<ret[i].length;j++){
				ret[i][j]=getZ(png.getRGB(i, j));
			}
		}
		return ret;
	}

	/**
	 * 標高の配列（double[][]）を標高PNGに変換
	 * @param z
	 * @return
	 */
	public static BufferedImage doubleToImg(double[][] z){
		BufferedImage ret=new BufferedImage(z.length,z[0].length,BufferedImage.TYPE_INT_RGB);
		for(int i=0;i<z.length;i++){
			for(int j=0;j<z[i].length;j++){
				ret.setRGB(i, j, getRGB(z[i][j]));
			}
		}
		return ret;
	}

	/**
	 * 標高PNGの減算
	 * @param b1 標高PNG1
	 * @param b2 標高PNG2
	 * @return 標高PNG（標高PNG1-標高PNG2）
	 * @throws IllegalArgumentException
	 */
	public static BufferedImage subZImage(BufferedImage b1,BufferedImage b2)throws IllegalArgumentException{
		int w=b1.getWidth();
		int h=b1.getHeight();
		if(w==b2.getWidth()&&h==b2.getHeight()){
			BufferedImage ret=new BufferedImage(w,h,BufferedImage.TYPE_INT_RGB);
			for(int i=0;i<w;i++){
				for(int j=0;j<h;j++){
					double v1=getZ(b1.getRGB(i, j));
					double v2=getZ(b2.getRGB(i, j));
					if(Double.isNaN(v1)||Double.isNaN(v2)){
						ret.setRGB(i, j, NA);
					}
					double v3=v1-v2;
					if(v3<0){
						ret.setRGB(i, j, NA);
					}else{
						ret.setRGB(i, j, getRGB(v3));
					}
				}
			}
			return ret;
		}else{
			throw new IllegalArgumentException("Image size is different");
		}
	}

	/**
	 * 標高PNGデータをテキストに出力
	 * @param png 標高PNG
	 * @param af アフィン変換
	 * @param out 出力ファイル
	 * @throws IOException
	 */
	public static void output(BufferedImage png,AffineTransform af,File out) throws IOException{
		Point2D dst=new Point2D.Double();
		BufferedWriter bw=new BufferedWriter(new FileWriter(out));
		int n=1;
		int w=png.getWidth();
		int h=png.getHeight();
		for(int i=0;i<w;i++){
			for(int j=0;j<h;j++){
				dst=af.transform(new Point2D.Double(i,j), dst);
				double v=getZ(png.getRGB(i, j));
				if(Double.isNaN(h))continue;
				bw.write(Integer.toString(n++)+","+dst.getX()+","+dst.getY()+","+Double.toString(v)+"\n");
			}
			bw.flush();
		}
		bw.close();
	}
}
~~~

## Lasデータセットの読み込み
Lasデータセットの読み込みには「[laszip4j](https://github.com/mreutegg/laszip4j)」を使用しました。
「[laszip4j](https://github.com/mreutegg/laszip4j)」のライセンスは「GNU Lesser General Public License v2.1」です。

複数のLasファイルを読み込んで結合し、点群から航空写真とDSMを生成するLASFileReaderクラスを作成しました。
なお、これら点群データは平面直角第8系の座標系のようです。

~~~java:LASFileReader
import java.awt.Color;
import java.awt.geom.AffineTransform;
import java.awt.geom.Point2D;
import java.awt.geom.Rectangle2D;
import java.awt.image.BufferedImage;
import java.io.File;
import java.io.IOException;
import java.util.ArrayList;
import java.util.Calendar;
import java.util.Date;
import java.util.List;

import javax.imageio.ImageIO;

import com.github.mreutegg.laszip4j.LASHeader;
import com.github.mreutegg.laszip4j.LASPoint;
import com.github.mreutegg.laszip4j.LASReader;

import net.termat.geo.ElevationPngUtil;

public class LASFileReader {
	private LASReader reader;
	private Rectangle2D bounds;
	private AffineTransform af;
	private BufferedImage img;
	private BufferedImage dem;
	private double[] minmaxZ;
	private Date date;
	private double xScalefactor;
	private double xOffset;
	private double yScalefactor;
	private double yOffset;
	private double zScalefactor;
	private double zOffset;
	private float basef=65536f;

	/**
	 *
	 * @param f Lasファイルの配列
	 * @param mPerPixel 1ピクセルの長さ(m)
	 * @throws Exception
	 */
	public static void outputLasDataUnion(File[] f,double mPerPixel)throws Exception{
		List<LASFileReader> list=new ArrayList<LASFileReader>();
		Rectangle2D rect=null;
		for(int i=0;i<f.length;i++){
			if(i==0){
				LASFileReader la=new LASFileReader(f[i],mPerPixel);
				rect=la.getBounds();
				list.add(la);
			}else{
				LASFileReader la=new LASFileReader(f[i],mPerPixel);
				rect=rect.createUnion(la.getBounds());
				list.add(la);
			}
		}
		AffineTransform af=new AffineTransform(new double[]{mPerPixel,0,0,-mPerPixel,rect.getX(),rect.getY()+rect.getHeight()});
		double width=Math.ceil(rect.getWidth())/mPerPixel;
		double height=Math.ceil(rect.getHeight())/mPerPixel;
		BufferedImage img=new BufferedImage((int)width,(int)height,BufferedImage.TYPE_INT_RGB);
		BufferedImage dem=new BufferedImage((int)width,(int)height,BufferedImage.TYPE_INT_RGB);
		for(int i=0;i<dem.getWidth();i++){
			for(int j=0;j<dem.getHeight();j++){
				dem.setRGB(i, j, ElevationPngUtil.NA);
			}
		}
		for(LASFileReader al : list){
			al.createImage(img,dem,af.createInverse());
		}
		String path=f[0].getAbsolutePath();
		File pf=new File(path.replace(".las", "_photp.jpg"));
		File df=new File(path.replace(".las", "_org.png"));
		ImageIO.write(img, "jpg", pf);
		ImageIO.write(dem, "png", df);
	}

	/**
	 *
	 * @param f Lasファイル
	 * @param mPerPixel 1ピクセルの長さ(m)
	 * @throws Exception
	 */
	public LASFileReader(File f,double mPerPixel)throws IOException{
		reader = new LASReader(f);
		readHeader();
		af=new AffineTransform(new double[]{mPerPixel,0,0,-mPerPixel,bounds.getX(),bounds.getY()+bounds.getHeight()});
		double width=Math.ceil(bounds.getWidth())/mPerPixel;
		double height=Math.ceil(bounds.getHeight())/mPerPixel;
		img=new BufferedImage((int)width,(int)height,BufferedImage.TYPE_INT_RGB);
		dem=new BufferedImage((int)width,(int)height,BufferedImage.TYPE_INT_RGB);
		for(int i=0;i<(int)width;i++){
			for(int j=0;j<(int)height;j++){
				dem.setRGB(i, j, ElevationPngUtil.NA);
			}
		}
	}

	/**
	 * Lasデータ領域の取得
	 * @return
	 */
	public Rectangle2D getBounds() {
		return bounds;
	}

	/*
	 * ヘッダーの読み込み
	 */
	private void readHeader(){
		LASHeader h=reader.getHeader();
		bounds=new Rectangle2D.Double(h.getMinX(),h.getMinY(),h.getMaxX()-h.getMinX(),h.getMaxY()-h.getMinY());
		minmaxZ=new double[]{h.getMinZ(),h.getMaxZ()};
		Calendar cal=Calendar.getInstance();
		cal.set(Calendar.YEAR, (int)h.getFileCreationYear());
		cal.set(Calendar.DAY_OF_YEAR, (int)h.getFileCreationDayOfYear());
		date=cal.getTime();
		xOffset=h.getXOffset();
		xScalefactor=h.getXScaleFactor();
		yOffset=h.getYOffset();
		yScalefactor=h.getYScaleFactor();
		zOffset=h.getZOffset();
		zScalefactor=h.getZScaleFactor();
	}

	/**
	 *  画像の生成
	 * @param img 航空写真
	 * @param dem 標高PNG画像
	 * @param at  アフィン変換
	 */
	public void createImage(BufferedImage img,BufferedImage dem,AffineTransform at){
	    for (LASPoint p : reader.getPoints()) {
	    	double x=xScalefactor*(double)p.getX()+xOffset;
	    	double y=yScalefactor*(double)p.getY()+yOffset;
	    	double z=zScalefactor*(double)p.getZ()+zOffset;
	    	Point2D px=at.transform(new Point2D.Double(x,y), new Point2D.Double());
	    	float r=((float)p.getRed())/basef;
	    	float g=((float)p.getGreen())/basef;
	    	float b=((float)p.getBlue())/basef;
	    	int col=new Color(r,g,b).getRGB();
	    	img.setRGB((int)px.getX(), (int)px.getY(), col);
	    	dem.setRGB((int)px.getX(), (int)px.getY(), ElevationPngUtil.getRGB(z));
	    }
	}

	/**
	 * 画像の生成
	 */
	public void createImage(){
		AffineTransform at=null;
		try{
			at=af.createInverse();
		}catch(Exception e){e.printStackTrace();}
	    for (LASPoint p : reader.getPoints()) {
	    	double x=xScalefactor*(double)p.getX()+xOffset;
	    	double y=yScalefactor*(double)p.getY()+yOffset;
	    	double z=zScalefactor*(double)p.getZ()+zOffset;
	    	Point2D px=at.transform(new Point2D.Double(x,y), new Point2D.Double());
	    	float r=((float)p.getRed())/basef;
	    	float g=((float)p.getGreen())/basef;
	    	float b=((float)p.getBlue())/basef;
	    	int col=new Color(r,g,b).getRGB();
	    	img.setRGB((int)px.getX(), (int)px.getY(), col);
	    	dem.setRGB((int)px.getX(), (int)px.getY(), ElevationPngUtil.getRGB(z));
	    }
	}
}
~~~

## ■結果
[静岡県ポイントクラウドDB](https://pointcloud.pref.shizuoka.jp/)よりダウンロードした点群データ（今回は30XXX00010025-1.las～30XXX00010025-6.lasのデータを使用）をLASFileReaderクラスで処理すると、以下の航空写真と標高PNG画像が生成されます。
![01.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/55735/a1428bd3-b72a-010b-aad7-6dae751cad41.jpeg)

## ■今後
[静岡県ポイントクラウドDB](https://pointcloud.pref.shizuoka.jp/)の様な公開データがあると、市民レベルでも色々な取り組みに活用できるように思います。

なお、今回、山林部の点群データを使用したのは、国土地理院で紹介されている「[航空レーザ測量データを用いた樹高等のデータ作成](https://www.gsi.go.jp/chirijoho/chirijoho40069.html)」を参考に森林部の樹冠解析に挑戦してみたいと考えたからです。

→[Javaで静岡県ポイントクラウドDBのデータを読み込んで、樹高等を検出してみる。](https://qiita.com/t-mat/items/f50c255e80c70705a702)