Read the data of Shizuoka point cloud DB with Java and generate aerial photograph and elevation PNG.
■ Overview
Recently, I am studying geographic information processing with half my work and half my hobbies. When I was searching for open point cloud data, I came to Shizuoka Point Cloud DB. In Shizuoka Point Cloud DB, the data is published in the Las dataset, so read Las in Java and aerial photographs and elevation data images (elevation PNG). ) Was generated.
■ Altitude PNG processing
Elevation PNG used in the Geographical Survey Institute elevation tile was used to convert the elevation into data. The paper at altitude PNG is as follows. → "[PNG Elevation Tile-Consideration and Implementation of Elevation File Format Suitable for Web Use- (Yoshiharu Nishioka / Juri Nagatsu, Information Geoinformatics Vol. 26, No. 4, 2015)](https://www.jstage.jst. go.jp/article/geoinformatics/26/4/26_155/_pdf) "
The ElevationPngUtil class was created with reference to the Geographical Survey Institute website and thesis.
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;
/**
* Elevation PNG utility class
*/
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(){}
/ * Elevation PNG NA value * /
public static int NA=P23;
/ * Convert RGB values to Int * /
private static int rgb2Int(int[] c){
return new Color(c[0],c[1],c[2]).getRGB();
}
/**
* Convert elevation value to elevation PNG pixel value
* @param z Elevation value (m)
* @return Elevation PNG pixel value
*/
public static int getRGB(double z){
if(Double.isNaN(z))return P23;
return rgb2Int(getRGBColor(z));
}
/ * Convert elevation value to 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};
}
/**
* Convert elevation PNG pixel value to elevation value
* @param intColor Elevation PNG pixel value
* @return Altitude value (m)
*/
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};
}
/**
* Convert elevation PNG to an array of elevations (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;
}
/**
* Convert an array of elevations (double [] []) to elevation 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;
}
/**
* Subtraction of elevation PNG
* @ param b1 Altitude PNG1
* @ param b2 Altitude PNG2
* @return Altitude PNG (elevation PNG1-elevation 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");
}
}
/**
* Output elevation PNG data to text
* @param png Altitude PNG
* @param af Affine transformation
* @param out Output file
* @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();
}
}
Loading Las dataset
I used "laszip4j" to load the Las dataset. The license for "laszip4j" is "GNU Lesser General Public License v2.1".
I created a LASFileReader class that reads and combines multiple Las files and generates aerial photographs and DSMs from point clouds. In addition, these point cloud data seems to be the coordinate system of the plane right angle 8th system.
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 file array
* @param mPerPixel 1 pixel length (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 file
* @param mPerPixel 1 pixel length (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);
}
}
}
/**
* Get Las data area
* @return
*/
public Rectangle2D getBounds() {
return bounds;
}
/*
* Read header
*/
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();
}
/**
* Image generation
* @param img Aerial photography
* @param dem elevation PNG image
* @param at Affine transformation
*/
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));
}
}
/**
* Image generation
*/
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));
}
}
}
■ Results
Point cloud data downloaded from Shizuoka Point Cloud DB (this time, data from 30XXX00010025-1.las to 30XXX00010025-6.las is used) with LASFileReader class. The process will generate the following aerial and elevation PNG images.
■ Future
If there is public data such as Shizuoka Point Cloud DB, I think that it can be used for various efforts even at the citizen level.
This time, the point cloud data of the forest area was used in "Creation of data such as tree height using aerial laser survey data" introduced by the Geospatial Information Authority of Japan. This is because I wanted to try the canopy analysis of the forest area with reference to "jp / chirijoho / chirijoho40069.html)".
→ Read the data of Shizuoka point cloud DB with Java and try to detect the tree height etc.