■概要
最近、仕事半分、趣味半分で地理情報処理の勉強をしています。
その中でオープンな点群データがないか検索していたら、静岡県ポイントクラウドDBに行きつきました。
静岡県ポイントクラウドDBでは、Lasデータセットでデータが公開されているので、JavaでLasを読み込んで航空写真と標高データ画像（標高PNG）を生成してみました。
■標高PNG処理
標高のデータ化には国土地理院標高タイルで使用されている標高PNGを用いました。
標高PNGの論文はこちらの様です。
→「PNG標高タイル －Web利用に適した標高ファイルフォーマットの考察と実装－（西岡芳晴・長津樹理、情報地質第26巻第4号、2015）」
地理院HP及び論文を参考に、ElevationPngUtilクラスを作成しました。
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」を使用しました。
「laszip4j」のライセンスは「GNU Lesser General Public License v2.1」です。
複数のLasファイルを読み込んで結合し、点群から航空写真とDSMを生成するLASFileReaderクラスを作成しました。
なお、これら点群データは平面直角第8系の座標系のようです。
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よりダウンロードした点群データ（今回は30XXX00010025-1.las～30XXX00010025-6.lasのデータを使用）をLASFileReaderクラスで処理すると、以下の航空写真と標高PNG画像が生成されます。
■今後
静岡県ポイントクラウドDBの様な公開データがあると、市民レベルでも色々な取り組みに活用できるように思います。
なお、今回、山林部の点群データを使用したのは、国土地理院で紹介されている「航空レーザ測量データを用いた樹高等のデータ作成」を参考に森林部の樹冠解析に挑戦してみたいと考えたからです。