1.概要
先週末、兵庫県が1mメッシュの全県域3次元データセットを公開をしました。
公開データはDSM、DEM、CS立体図で、ライセンスはC.C.4.0、「目的を問わず二次利用可能なデータですので、様々な用途でご利用ください。」とのことです。
これまで公開されていた三次元データセットとしては、国土地理院：基盤地図情報サイトの数値標高モデル（5ｍメッシュDEM）がありましたが、兵庫県全域数値地形図はその25倍のデータ密度であり、DEMに加え、建物や森林を含めたDSMが公開されています。
さっそくこれを利用して、都市や山林部の三次元モデルを作りたいと思いました。
2.今回作るもの
兵庫県公開のDSMは、平面直角座標第5系のx,y座標と標高値からなるデータセットです。三次元モデルを作るには色情報が欲しいところです。
このため、まず、地理院タイル：全国最新写真の画像タイルから各メッシュの色情報を取得し、DSMをXYZRGB形式のデータに変換するツールを作成します。
次いで、XYZRGBの頂点データから面を生成し、PLY形式の三次元モデルを生成するツールを作成します。
なお、PLY形式の仕様については、下記のサイトを参考にしました。
3.DSM ->XYZRGB変換ツールの実装
(1)処理の流れ
処理の流れは、以下のとおりです。
1.DSMのx,y座標(平面直角座標第5系）を検査し、DSM点群の領域を取得する。
2.DSM点群領域を緯度経度の矩形領域に座標変換する。
3.緯度経度の矩形領域をピクセル座標に変換し、取得するタイル座標を把握する。
4.DSMに対応した座標の画像タイルを取得し、連結する。
5.DSM点群座標をピクセル座標に変換し、取得した航空写真画像から対応するピクセルの色情報を取得する。
6.「x y z r g b」の様式でテキストファイルに出力する。
(2)座標変換クラス（平面直角座標系<->緯度経度）の実装
まず、平面直角座標系のxy座標を緯度経度に変換するクラスを実装します。
国土地理院時報の「Gauss-Krüger 投影における経緯度座標及び平面直角座標相互間の座標換算についてのより簡明な計算方法」を参考に下記のクラスを作成しました。
LonlatXYT.java
import java.awt.geom.Point2D;
public class LonLatXY {
private static final double a=6378137;
private static final double rf=298.257222101;
private static final double m0=0.9999;
private static final double s2r=Math.PI/648000;
private static final double n=0.5/(rf-0.5);
private static final double n15=1.5*n;
private static final double anh=0.5*a/(1+n);
private static final double nsq=n*n;
private static final double e2n=2*Math.sqrt(n)/(1+n);
private static final double ra=2*anh*m0*(1+nsq/4+nsq*nsq/64);
private static int jt=5;
private static int jt2=2*jt;
private static double ep=1.0;
private static double[] e=getE();
private static final double[] phi0=new double[]{0,33,33,36,33,36,36,36,36,36,40,44,44,44,26,26,26,26,20,26};
private static final double[] lmbd0=new double[]{0,7770,7860,7930,8010,8060,8160,8230,8310,8390,8450,8415,8535,8655,8520,7650,7440,7860,8160,9240};
private static double[] alp=getAlp();
private static double[] beta=getBeta();
private static double[] dlt=getDlt();
private static double[] getAlp(){
double[] alp=new double[6];
alp[1]=(1.0/2.0+(-2.0/3.0+(5.0/16.0+(41.0/180.0-127.0/288.0*n)*n)*n)*n)*n;
alp[2]=(13.0/48.0+(-3.0/5.0+(557.0/1440.0+281.0/630.0*n)*n)*n)*nsq;
alp[3]=(61.0/240.0+(-103.0/140.0+15061.0/26880.0*n)*n)*n*nsq;
alp[4]=(49561.0/161280.0-179.0/168.0*n)*nsq*nsq;
alp[5]=34729.0/80640.0*n*nsq*nsq;
return alp;
}
private static double[] getBeta(){
double[] beta=new double[6];
beta[1]=(1.0/2.0+(-2.0/3.0+(37.0/96.0+(-1.0/360.0-81.0/512.0*n)*n)*n)*n)*n;
beta[2]=(1.0/48.0+(1.0/15.0+(-437.0/1440.0+46.0/105.0*n)*n)*n)*nsq;
beta[3]=(17.0/480.0+(-37.0/840.0-209.0/4480.0*n)*n)*n*nsq;
beta[4]=(4397.0/161280.0-11.0/504.0*n)*nsq*nsq;
beta[5]=4583.0/161280.0*n*nsq*nsq;
return beta;
}
private static double[] getDlt(){
double[] dlt=new double[7];
dlt[1]=(2.0+(-2.0/3.0+(-2.0+(116.0/45.0+(26.0/45.0-2854.0/675.0*n)*n)*n)*n)*n)*n;
dlt[2]=(7.0/3.0+(-8.0/5.0+(-227.0/45.0+(2704.0/315.0+2323.0/945.0*n)*n)*n)*n)*nsq;
dlt[3]=(56.0/15.0+(-136.0/35.0+(-1262.0/105.0+73814.0/2835.0*n)*n)*n)*n*nsq;
dlt[4]=(4279.0/630.0+(-332.0/35.0-399572.0/14175.0*n)*n)*nsq*nsq;
dlt[5]=(4174.0/315.0-144838.0/6237.0*n)*n*nsq*nsq;
dlt[6]=601676.0/22275.0*nsq*nsq*nsq;
return dlt;
}
private static double[] getE(){
double[] e=new double[jt2+2];
for(int k=1;k<=jt;k++){
ep*=e[k]=n15/k-n;
e[k+jt]=n15/(k+jt)-n;
}
return e;
}
public static Point2D xyToLonLat(int num,double xx,double yy){
double x=yy;
double y=xx;
double xi=(x+m0*Merid(2*phi0[num]*3600*s2r))/ra;
double xip=xi;
double eta=y/ra;
double etap=eta;
double sgmp=1;
double taup=0;
for(int j=beta.length-1;j>0;j--){
double besin=beta[j]*Math.sin(2*j*xi);
double becos=beta[j]*Math.cos(2*j*xi);
xip -=besin*Math.cosh(2*j*eta);
etap -=becos*Math.sinh(2*j*eta);
sgmp -=2*j*becos*Math.cosh(2*j*eta);
taup +=2*j*besin*Math.sinh(2*j*eta);
}
double sxip=Math.sin(xip);
double cxip=Math.cos(xip);
double shetap=Math.sinh(etap);
double chetap=Math.cosh(etap);
double chi=Math.asin(sxip/chetap);
double phi=chi;
for(int j=dlt.length-1;j>=0;j--){
phi +=dlt[j]*Math.sin(2*j*chi);
}
double nphi=(1-n)/(1+n)*Math.tan(phi);
double lmbd=lmbd0[num]*60+Math.atan2(shetap, cxip)/s2r;
double lat=phi/s2r/3600;
double lon=lmbd/3600;
return new Point2D.Double(lon,lat);
}
public static Point2D lonlatToXY(int num,double lon,double lat){
double phirad=Math.toRadians(lat);
double lmbddeg=Math.floor(lon);
double lmbdmin=Math.floor(60.0*(lon-lmbddeg));
double lmbdsec=lmbddeg*3600.0+lmbdmin*60.0+(lon-lmbddeg-lmbdmin/60)*3600.0;
double sphi=Math.sin(phirad);
double nphi=(1-n)/(1+n)*Math.tan(phirad);
double dlmbd=(lmbdsec-lmbd0[num]*60.0)*s2r;
double sdlmbd=Math.sin(dlmbd);
double cdlmbd=Math.cos(dlmbd);
double tchi=Math.sinh(atanh(sphi)-e2n*atanh(e2n*sphi));
double cchi=Math.sqrt(1+tchi*tchi);
double xip=Math.atan2(tchi, cdlmbd);
double xi=xip;
double etap=atanh(sdlmbd/cchi);
double eta=etap;
double sgm=1;
double tau=0;
for(int j=alp.length-1;j>0;j--){
double alsin=alp[j]*Math.sin(2*j*xip);
double alcos=alp[j]*Math.cos(2*j*xip);
xi +=alsin*Math.cosh(2*j*etap);
eta +=alcos*Math.sinh(2*j*etap);
sgm +=2*j*alcos*Math.cosh(2*j*etap);
tau +=2*j*alsin*Math.sinh(2*j*etap);
}
double x=ra*xi-m0*Merid(2*phi0[num]*3600*s2r);
double y=ra*eta;
return new Point2D.Double(x,y);
}
private static double Merid(double phi2) {
double dc=2.0*Math.cos(phi2);
double[] s=new double[jt2+2];
double[] t=new double[jt2+2];
s[1]=Math.sin(phi2);
for(int i=1;i<=jt2;i++){
s[i+1]=dc*s[i]-s[i-1];
t[i]=(1.0/i-4.0*i)*s[i];
}
double sum=0.0;
double c1=ep;
int j=jt;
while(j>0){
double c2=phi2;
double c3=2.0;
int l=j;
int m=0;
while(l>0){
c2 +=(c3/=e[l--])*t[++m]+(c3*=e[2*j-l])*t[++m];
}
sum +=c1*c1*c2 ; c1/=e[j--];
}
return anh*(sum+phi2);
}
private static double atanh(double v){
return 0.5*Math.log((1.0+v)/(1.0-v));
}
}
(3)航空写真タイル画像取得クラスの実装
地理院地図からタイル画像を取得・連結する下記クラスを作成しました。
GSITileReader.java
import java.awt.Graphics2D;
import java.awt.Shape;
import java.awt.geom.AffineTransform;
import java.awt.geom.GeneralPath;
import java.awt.geom.Point2D;
import java.awt.geom.Rectangle2D;
import java.awt.image.BufferedImage;
import java.io.BufferedWriter;
import java.io.File;
import java.io.FileWriter;
import java.io.IOException;
import java.net.URL;
import java.util.Arrays;
import java.util.HashSet;
import java.util.Set;
import javax.imageio.ImageIO;
public class GSITileReader {
private static final String base_url="https://cyberjapandata.gsi.go.jp/xyz/seamlessphoto/";
private static final double L=85.05112877980659;
public static BufferedImage getGSIImage(Rectangle2D plane,int map_num,int zoom,double scale)throws IOException{
Shape sp=getLonLatShapeAtPlane(map_num, plane.getX(), plane.getY(), plane.getWidth(), plane.getHeight());
Rectangle2D rect=sp.getBounds2D();
long[] topLeft=lonlatToPixel(zoom,rect.getX(),rect.getY());
long[] bottomRight=lonlatToPixel(zoom,rect.getX()+rect.getWidth(),rect.getY()+rect.getHeight());
Set<Long> x=new HashSet<Long>();
Set<Long> y=new HashSet<Long>();
for(long i=topLeft[1];i<=bottomRight[1];i++){
x.add((long)Math.ceil(i/256));
}
for(long i=bottomRight[2];i<=topLeft[2];i++){
y.add((long)Math.ceil(i/256));
}
Long[] xx=x.toArray(new Long[x.size()]);
Long[] yy=y.toArray(new Long[y.size()]);
Arrays.sort(xx);
Arrays.sort(yy);
BufferedImage im=new BufferedImage(xx.length*256,yy.length*256,BufferedImage.TYPE_INT_RGB);
Graphics2D g=im.createGraphics();
for(int i=0;i<xx.length;i++){
for(int j=0;j<yy.length;j++){
try{
String url=base_url+Integer.toString(zoom)+"/"+Long.toString(xx[i])+"/"+Long.toString(yy[j])+".jpg";
BufferedImage tmp=ImageIO.read(new URL(url));
g.drawImage(tmp, i*256, j*256, null);
}catch(IOException e){
e.printStackTrace();
}
}
}
g.dispose();
im=subImage(im,plane,(long)xx[0],(long)yy[0],map_num,zoom,scale);
return im;
}
private static BufferedImage subImage(BufferedImage src,Rectangle2D rect,long minX,long minY,int num,int zoom,double scale){
double xx=rect.getX();
double yy=rect.getY();
double ww=Math.abs(rect.getWidth());
double hh=Math.abs(rect.getHeight());
minX=minX*256;
minY=minY*256;
BufferedImage dst=new BufferedImage((int)(ww/scale),(int)(hh/scale),BufferedImage.TYPE_INT_RGB);
System.out.println(dst.getWidth()+"/"+dst.getHeight());
System.out.println(minX+"/"+minY);
for(int i=0;i<dst.getWidth();i++){
for(int j=0;j<dst.getHeight();j++){
double x=xx+scale*i;
double y=yy-scale*j;
Point2D p=LonLatXY.xyToLonLat(num, x, y);
long[] pc=lonlatToPixel(zoom, p.getX(), p.getY());
int px=(int)(pc[1]-minX);
int py=(int)(pc[2]-minY);
int color=src.getRGB(px, py);
dst.setRGB(i, j, color);
}
}
return dst;
}
public static long[] lonlatToPixel(int zoom,double lon,double lat){
long x=(long)(Math.pow(2, zoom+7)*(lon/180.0+1.0));
long y=(long)((Math.pow(2, zoom+7)/Math.PI)*(-atanh(Math.sin(Math.toRadians(lat)))+atanh(Math.sin(Math.toRadians(L)))));
return new long[]{(long)zoom,x,y};
}
private static double atanh(double v){
return 0.5*Math.log((1.0+v)/(1.0-v));
}
public static AffineTransform createTfwTransform(Rectangle2D rectXY,BufferedImage img){
double sx=rectXY.getWidth()/img.getWidth();
double sy=rectXY.getHeight()/img.getHeight();
double x=rectXY.getX();
double y=rectXY.getY()+rectXY.getHeight();
AffineTransform af=new AffineTransform(new double[]{sx,0,0,-sy,x,y});
return af;
}
public static Shape getLonLatShapeAtPlane(int num,double x,double y,double w,double h){
Point2D p1=LonLatXY.xyToLonLat(num, x, y);
Point2D p2=LonLatXY.xyToLonLat(num, x+w, y);
Point2D p3=LonLatXY.xyToLonLat(num, x+w, y+h);
Point2D p4=LonLatXY.xyToLonLat(num, x, y+h);
GeneralPath gp=new GeneralPath();
gp.moveTo(p1.getX(),p1.getY());
gp.lineTo(p2.getX(), p2.getY());
gp.lineTo(p3.getX(), p3.getY());
gp.lineTo(p4.getX(), p4.getY());
gp.closePath();
return gp;
}
public static void outTfw(AffineTransform af,File out)throws IOException{
BufferedWriter bw=new BufferedWriter(new FileWriter(out));
bw.write(af.getScaleX()+"\n");
bw.write(af.getShearX()+"\n");
bw.write(af.getShearY()+"\n");
bw.write(af.getScaleY()+"\n");
bw.write(af.getTranslateX()+"\n");
bw.write(af.getTranslateY()+"\n");
bw.close();
}
}
(4)XYZRGB出力クラス
最後にファイル入出力など、DSMをXYZRGBに変換するツールを作成しました。
本クラスでは、取得した航空写真画像からDSM座標に対応したピクセルのRGBを出力し、XYZRGB形式で出力しています。
一応、平面直角座標第5系以外にも対応できるようにしています。
XYZRgbCreator.java
import java.awt.BorderLayout;
import java.awt.Color;
import java.awt.Font;
import java.awt.datatransfer.DataFlavor;
import java.awt.datatransfer.Transferable;
import java.awt.datatransfer.UnsupportedFlavorException;
import java.awt.event.WindowAdapter;
import java.awt.event.WindowEvent;
import java.awt.geom.AffineTransform;
import java.awt.geom.Point2D;
import java.awt.geom.Rectangle2D;
import java.awt.image.BufferedImage;
import java.io.BufferedReader;
import java.io.BufferedWriter;
import java.io.File;
import java.io.FileReader;
import java.io.FileWriter;
import java.io.IOException;
import java.util.List;
import javax.imageio.ImageIO;
import javax.swing.JComboBox;
import javax.swing.JFrame;
import javax.swing.JLabel;
import javax.swing.JOptionPane;
import javax.swing.JPanel;
import javax.swing.SwingUtilities;
import javax.swing.TransferHandler;
import javax.swing.UIManager;
public class XYZRgbCreator {
private JFrame frame;
private JComboBox<String> zone;
public XYZRgbCreator(){
frame=new JFrame();
frame.setTitle("XYZRGB");
frame.setDefaultCloseOperation(JFrame.DO_NOTHING_ON_CLOSE);
try {
UIManager.setLookAndFeel("com.sun.java.swing.plaf.nimbus.NimbusLookAndFeel");
SwingUtilities.updateComponentTreeUI(frame);
}catch(Exception e){
try {
UIManager.setLookAndFeel("com.sun.java.swing.plaf.windows.WindowsLookAndFeel");
SwingUtilities.updateComponentTreeUI(frame);
}catch(Exception ee){
ee.printStackTrace();
}
}
WindowAdapter wa=new WindowAdapter(){
@Override
public void windowClosing(WindowEvent e) {
close();
}
};
frame.addWindowListener(wa);
frame.getContentPane().setLayout(new BorderLayout());
zone=new JComboBox<String>(new String[]{
null,"平面直角第01系","平面直角第02系","平面直角第03系","平面直角第04系","平面直角第05系","平面直角第06系","平面直角第07系","平面直角第08系",
"平面直角第09系","平面直角第10系","平面直角第11系","平面直角第12系","平面直角第13系","平面直角第14系","平面直角第15系","平面直角第16系",
"平面直角第17系","平面直角第18系","平面直角第19系"
});
zone.setSelectedIndex(5);
frame.getContentPane().add(zone,BorderLayout.NORTH);
JLabel label=new JLabel("DSM > XYZRGB");
label.setFont(new Font(Font.SANS_SERIF,Font.BOLD,24));
label.setVerticalAlignment(JLabel.CENTER);
label.setHorizontalAlignment(JLabel.CENTER);
JPanel jp=new JPanel(new BorderLayout());
jp.add(label,BorderLayout.CENTER);
frame.getContentPane().add(jp,BorderLayout.CENTER);
frame.setSize(480,480);
frame.setResizable(false);
DropFileHandler handler=new DropFileHandler();
label.setTransferHandler(handler);
jp.setTransferHandler(handler);
}
private void close(){
int id=JOptionPane.showConfirmDialog(frame, "Exit?", "Info", JOptionPane.YES_NO_OPTION,JOptionPane.INFORMATION_MESSAGE);
if(id==JOptionPane.YES_OPTION){
frame.setVisible(false);
System.exit(0);
}
}
private class DropFileHandler extends TransferHandler {
private static final long serialVersionUID = 1L;
@Override
public boolean canImport(TransferSupport support) {
if (!support.isDrop()) {
return false;
}
if (!support.isDataFlavorSupported(DataFlavor.javaFileListFlavor)) {
return false;
}
return true;
}
@SuppressWarnings("unchecked")
@Override
public boolean importData(TransferSupport support) {
if (!canImport(support)) {
return false;
}
Transferable t = support.getTransferable();
try {
List<File> files = (List<File>) t.getTransferData(DataFlavor.javaFileListFlavor);
for (File file : files){
if(file.isDirectory())continue;
try{
process(file);
}catch(IOException e){
e.printStackTrace();
}
}
} catch (UnsupportedFlavorException | IOException e) {
e.printStackTrace();
}
return true;
}
}
private void process(File f) throws IOException{
if(!(f.getName().endsWith(".txt")||f.getName().endsWith(".xyz")))return;
frame.repaint();
int bwz=zone.getSelectedIndex();
Rectangle2D rectXY=getBounds(f);
BufferedImage img=GSITileReader.getGSIImage(rectXY,bwz,18,0.5);
String name=f.getName().substring(0,f.getName().lastIndexOf("."));
File dir=f.getParentFile();
ImageIO.write(img, "jpg", new File(dir.getAbsolutePath()+"/"+name+".jpg"));
AffineTransform af=new AffineTransform(new double[]{
0.5,0,0,-0.5,rectXY.getX(),rectXY.getY()});
GSITileReader.outTfw(af, new File(dir.getAbsolutePath()+"/"+name+".jgw"));
try{
af=af.createInverse();
BufferedReader br=new BufferedReader(new FileReader(f));
File out=new File(dir.getAbsolutePath()+"/"+name+"_color.txt");
BufferedWriter bw=new BufferedWriter(new FileWriter(out));
String line=null;
String str=null;
Rectangle2D imgRect=new Rectangle2D.Double(0, 0, img.getWidth(),img.getHeight());
while((line=br.readLine())!=null){
line=line.replaceAll(" ",",");
String[] sp=line.split(",");
double x=Double.parseDouble(sp[0]);
double y=Double.parseDouble(sp[1]);
Point2D p=af.transform(new Point2D.Double(x,y), new Point2D.Double());
int xx=(int)Math.floor(p.getX());
int yy=(int)Math.floor(p.getY());
if(imgRect.contains(xx, yy)){
int col=img.getRGB(xx, yy);
Color color=new Color(col);
str=line+" "+Integer.toString(color.getRed())+" "+Integer.toString(color.getGreen())+" "+Integer.toString(color.getBlue())+"\n";
}else{
str=line+" 0 0 0\n";
}
bw.write(str);
bw.flush();
}
br.close();
bw.close();
}catch(Exception e){
e.printStackTrace();
}
}
private Rectangle2D getBounds(File f)throws IOException{
double xmin=Double.MAX_VALUE;
double xmax=-Double.MAX_VALUE;
double ymin=Double.MAX_VALUE;
double ymax=-Double.MAX_VALUE;
BufferedReader br=new BufferedReader(new FileReader(f));
String line=null;
while((line=br.readLine())!=null){
line=line.replaceAll(" ",",");
String[] sp=line.split(",");
double x=Double.parseDouble(sp[0]);
double y=Double.parseDouble(sp[1]);
xmin=Math.min(x, xmin);
xmax=Math.max(x, xmax);
ymin=Math.min(y, ymin);
ymax=Math.max(y, ymax);
}
br.close();
Rectangle2D ret=new Rectangle2D.Double(xmin,ymax,xmax-xmin+1,-(ymax-ymin)-1);
return ret;
}
public static void main(String[] args){
XYZRgbCreator gp=new XYZRgbCreator();
gp.frame.setLocationRelativeTo(null);
gp.frame.setVisible(true);
}
}
(5)成果物
XYZRgbCreatorを実行すると下記のウインドウが表示されます。兵庫県全域数値地形図DSMファイルをウインドウ内にドロップすれば、処理が開始され、xyzrgbファイルが出力されます。
なお、本来は処理の進捗情報を表示すべきなのですが、今回は割愛しました。
4.RGBXYZ -> PLY変換ツールの実装
次いで、XYZRGBファイルをPLYファイルに変換するツールを実装します。
DSMモデルはサーフェスモデルと呼ばれる表面のみが定義された3次元構造のモデルです。このため、2次元Delaunay分割を実行し、三角形ポリゴンをFaceにすれば三次元モデルが生成できると考えました。
なお、Delaunay分割には、2D Delaunay Triangulationライブラリ「TINFOUR」を利用しました。
mvn repositiry:tinfourから必要な情報を確認し、pomファイルのdependancisに追記しました。
(1)処理の流れ
処理の流れは以下のとおりです。
1.XYZRGBファイルのXYZ情報から頂点リストを、RGB情報から色情報リストを作成する。
2.IncrementalTinインスタンスに頂点リストを渡し、2次元Delaunay分割を実行する。
3.PLYの仕様に合わせ、必要事項を出力する。
4.頂点と色情報を出力する。
5.面（Face)情報を出力する。
(2)PlyCreatorクラスの実装
上記の処理を実行するPlyCreatorクラスを実装します。
なお、めんどくさくなったので、コマンドラインツールとして実装しています。
PlyCreator.java
import java.awt.Color;
import java.io.BufferedReader;
import java.io.BufferedWriter;
import java.io.File;
import java.io.FileReader;
import java.io.IOException;
import java.nio.charset.Charset;
import java.nio.file.Files;
import java.util.ArrayList;
import java.util.List;
import java.util.Observable;
import java.util.function.Consumer;
import org.tinfour.common.SimpleTriangle;
import org.tinfour.common.Vertex;
import org.tinfour.standard.IncrementalTin;
import org.tinfour.utils.TriangleCollector;
public class PlyCreator extends Observable{
private List<Vertex> vertexs;
private List<Color> colors;
private IncrementalTin tin;
public PlyCreator(){}
private void outVertex(BufferedWriter bw)throws IOException{
System.out.println("Vertex出力");
for(Vertex v : vertexs){
StringBuffer buf=new StringBuffer();
buf.append(Float.toString((float)v.getX())+" ");
buf.append(Float.toString((float)v.getY())+" ");
buf.append(Float.toString((float)v.getZ())+" ");
Color c=colors.get(v.getIndex());
buf.append(Integer.toString(c.getRed())+" ");
buf.append(Integer.toString(c.getGreen())+" ");
buf.append(Integer.toString(c.getBlue())+"\n");
writeBytes(bw,buf.toString());
}
}
private void outFace(BufferedWriter bw)throws IOException{
System.out.println("Face出力");
Consumer<SimpleTriangle> cons=new Consumer<SimpleTriangle>() {
@Override
public void accept(SimpleTriangle arg){
StringBuffer buf=new StringBuffer();
buf.append("3");
buf.append(" "+Integer.toString(arg.getVertexA().getIndex()));
buf.append(" "+Integer.toString(arg.getVertexB().getIndex()));
buf.append(" "+Integer.toString(arg.getVertexC().getIndex()));
buf.append("\n");
try{
writeBytes(bw,buf.toString());
}catch(IOException e){
e.printStackTrace();
}
}
};
TriangleCollector.visitSimpleTriangles(tin, cons);
}
private void writeBytes(BufferedWriter bw,String str)throws IOException{
bw.write(str,0,str.length());
}
public void writePLY(File f)throws IOException{
System.out.println("PLY出力");
Charset charset = Charset.forName("US-ASCII");
BufferedWriter bw = Files.newBufferedWriter(f.toPath(),charset);
writeBytes(bw,"ply\n");
writeBytes(bw,"format ascii 1.0\n");
writeBytes(bw,"element vertex "+Integer.toString(vertexs.size())+"\n");
writeBytes(bw,"property float x\n");
writeBytes(bw,"property float y\n");
writeBytes(bw,"property float z\n");
writeBytes(bw,"property uchar red\n");
writeBytes(bw,"property uchar green\n");
writeBytes(bw,"property uchar blue\n");
writeBytes(bw,"element face "+Integer.toString(tin.countTriangles().getCount())+"\n");
writeBytes(bw,"property list uchar int vertex_index\n");
writeBytes(bw,"end_header\n");
outVertex(bw);
outFace(bw);;
bw.close();
}
public void readXYZGRB(File f,String separator)throws IOException{
System.out.println("ファイル読み込み");
BufferedReader br=new BufferedReader(new FileReader(f));
String line=null;
vertexs=new ArrayList<Vertex>();
colors=new ArrayList<Color>();
int id=0;
while((line=br.readLine())!=null){
String[] sp=line.split(separator);
Vertex v=new Vertex(
Double.parseDouble(sp[0]),
Double.parseDouble(sp[1]),
Double.parseDouble(sp[2]),
id++
);
vertexs.add(v);
Color c=new Color(
Integer.parseInt(sp[3]),
Integer.parseInt(sp[4]),
Integer.parseInt(sp[5])
);
colors.add(c);
}
br.close();
System.out.println("TIN処理");
tin=new IncrementalTin();
tin.add(vertexs, null);
}
public static void main(String[] args){
PlyCreator pc=new PlyCreator();
File in=new File(args[0]);
File out=new File(args[1]);
try{
pc.readXYZGRB(in, ",");
pc.writePLY(out);
}catch(Exception e){
e.printStackTrace();
}
}
}
(3)成果物
入力ファイル(XYZRGBファイル）と出力ファイル（PLYファイル）を引数としてPlyCreatorクラスを実行するとPLY形式の三次元モデルが生成されます。
MeshLab等の3次元ビュアーにPLYファイルを読み込ませると、作成したモデルが表示されます。
下の画像はJR三ノ宮駅付近のDSMから生成したモデル画像です。
5.まとめ
国土地理院航空写真から色情報を取得して兵庫県全域数値地形図DSMをXYZRGBファイルに変換し、PLY形式のモデルを生成してみました。
都市域のモデルをみると多少荒く感じられるかもしれませんが、以下の画像は佐用町の山林部のDSMをモデル化したものですが、1ｍメッシュなら林相による樹冠の違いが表現された地形モデルを得ることができます。
兵庫県ホームページでは「このデータや他のデータと最先端技術を組み合わせて、民・産・学・官の協働による地域課題の解決に向けた取組を推進するため、利活用のアイデア・提案を募集します。」とありますが、兵庫県全域数値地形図は、兵庫県全域のDSM、DEMが公開されているので、色々な事に使えそうです。
兵庫県以外では、静岡県が「静岡ポイントクラウドデータベース」で点群データを公開していますが、他県でもこうした公共測量データ公開の動きがあると面白いと思います。