LoginSignup
3
0

More than 3 years have passed since last update.

@t-mat

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

Last updated at Posted at 2019-09-16

■概要

最近、仕事半分、趣味半分で地理情報処理の勉強をしています。
その中でオープンな点群データがないか検索していたら、静岡県ポイントクラウド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画像が生成されます。
01.jpg

■今後

静岡県ポイントクラウドDBの様な公開データがあると、市民レベルでも色々な取り組みに活用できるように思います。

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

→Javaで静岡県ポイントクラウドDBのデータを読み込んで、樹高等を検出してみる。

3
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
3
0