基盤地図情報からLOD1のCityGMLを作ってみる。

#1. 概要
[基盤地図情報](https://www.gsi.go.jp/kiban/)の建物ポリゴンから、[citygml4j](https://github.com/citygml4j/citygml4j)を使って、LOD1のCityGMLを作ってみました。
建物の高さの情報は、1mメッシュのDSMとDEMが公開されている[兵庫県高精度3次元データ](https://web.pref.hyogo.lg.jp/kk26/hyogo-geo.html)を使用しています。
なお、ソースコード一式を[こちら](https://github.com/termat/citygml_lod1)にアップしています。
![001.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/55735/60920cb4-503e-c8de-5892-c56984f51539.jpeg)

#2.建物情報の抽出
今回は、神戸市三宮周辺のCityGMLを作ってみました。
[基盤地図情報](https://www.gsi.go.jp/kiban/)から三宮を含むデータをダウンロードし、[QGIS](https://qgis.org/ja/site/forusers/download.html)で対象区域の建物をクリップして平面直角第5系でGeojsonに出力しました。
この時、ついでに面積(AREA)と周囲長（PERIMETER）を属性を追加しています。
![003.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/55735/f00e4681-c055-8fbc-c963-34c975ebdb97.jpeg)


#3.DSM・DEMの取得
[兵庫県高精度3次元データ](https://web.pref.hyogo.lg.jp/kk26/hyogo-geo.html)のDSM、DEMタイルから、処理区域の標高データ（[PNG標高画像](https://maps.gsi.go.jp/development/demtile.html)）を取得します。
実装したクラスは以下のとおりです。
・WebTile.java：任意領域の地図画像タイルを収集・結合して出力するクラス。
・GeojsonProcesser：Geojson処理用クラス
・TileApp.java：DropしたGeojson領域の兵庫県DSM、DEMを出力するGUIツール。

```java:WebTile.java
package citygml_lod1.components;

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.net.URL;
import java.security.KeyManagementException;
import java.security.NoSuchAlgorithmException;
import java.security.SecureRandom;
import java.security.cert.CertificateException;
import java.security.cert.X509Certificate;
import java.util.HashMap;
import java.util.Map;

import javax.imageio.ImageIO;
import javax.net.ssl.HostnameVerifier;
import javax.net.ssl.HttpsURLConnection;
import javax.net.ssl.SSLContext;
import javax.net.ssl.SSLSession;
import javax.net.ssl.X509TrustManager;


public class WebTile {
	private static final double L=85.05112877980659;
	private String url;
	private Map<Point2D,BufferedImage> tiles;
	private BufferedImage img;
	private AffineTransform af;
	private int zoom;
	private double resolution;
	
	public WebTile(String url,int zoom,double res){
		this.url=url;
		this.zoom=zoom;
		this.resolution=res;
		tiles=new HashMap<>();
	}
	
	public void create(int coordSys,Rectangle2D xy)throws IOException{
		int w=(int)Math.ceil(xy.getWidth()/resolution);
		int h=(int)Math.ceil(xy.getHeight()/resolution);
		if(w%2==1)w++;
		if(h%w==1)h++;
		double[] param=new double[]{
				resolution,0,0,-resolution,xy.getX(),xy.getY()+xy.getHeight()};
		af=new AffineTransform(param);
		img=new BufferedImage(w,h,BufferedImage.TYPE_INT_RGB);
		for(int i=0;i<w;i++){
			for(int j=0;j<h;j++){
				Point2D pxy=af.transform(new Point2D.Double(i,j),new Point2D.Double());
				Point2D lonlat=LonLatXY.xyToLonlat(coordSys, pxy.getX(), pxy.getY());
				Point2D pixel=lonlatToPixel(zoom,lonlat);
				Point2D tile=new Point2D.Double(Math.floor(pixel.getX()/256),Math.floor(pixel.getY()/256));
				if(tiles.containsKey(tile)){
					BufferedImage tmp=tiles.get(tile);
					if(tmp!=null){
						int xx=(int)pixel.getX()%256;
						int yy=(int)pixel.getY()%256;
						img.setRGB(i, j, tmp.getRGB(xx, yy));
					}
				}else{
					BufferedImage tmp=getTile(zoom,(long)tile.getX(),(long)tile.getY());
					if(tmp!=null){
						int xx=(int)pixel.getX()%256;
						int yy=(int)pixel.getY()%256;
						img.setRGB(i, j, tmp.getRGB(xx, yy));
					}
					tiles.put(tile, tmp);
				}
			}
		}
	}

	protected BufferedImage getTile(int zoom,long x,long y) {
		String uu=new String(url).replace("{z}", Integer.toString(zoom));
		uu=uu.replace("{x}", Long.toString(x));
		uu=uu.replace("{y}", Long.toString(y));
		try {
			HttpsURLConnection con=(HttpsURLConnection)new URL(uu).openConnection();
	        SSLContext sslContext = SSLContext.getInstance("SSL");
	        sslContext.init(null,
	                        new X509TrustManager[] { new LooseTrustManager() },
	                        new SecureRandom());

	        con.setSSLSocketFactory(sslContext.getSocketFactory());
	        con.setHostnameVerifier(new LooseHostnameVerifier());
	        BufferedImage tmp=ImageIO.read(con.getInputStream());
			if(tmp!=null)return tmp;
		}catch(IOException | KeyManagementException | NoSuchAlgorithmException e) {
			e.printStackTrace();
		}
		return null;
	}
	
	public BufferedImage getImage() {
		return img;
	}

	public AffineTransform getTransform() {
		return af;
	}

	private static Point2D lonlatToPixel(int zoom,Point2D p){
		long x=(long)(Math.pow(2, zoom+7)*(p.getX()/180.0+1.0));
		long y=(long)((Math.pow(2, zoom+7)/Math.PI)*(-atanh(Math.sin(Math.toRadians(p.getY())))+atanh(Math.sin(Math.toRadians(L)))));
		return new Point2D.Double(x,y);
	}

	private static double atanh(double v){
		return 0.5*Math.log((1.0+v)/(1.0-v));
	}
	
	private static class LooseTrustManager implements X509TrustManager {
	    @Override
	    public void checkClientTrusted(X509Certificate[] chain, String authType) throws CertificateException {}
	 
	    @Override
	    public void checkServerTrusted(X509Certificate[] chain, String authType) throws CertificateException {}
	 
	    @Override
	    public X509Certificate[] getAcceptedIssuers() {
	        return null;
	    }
	}
	
	private static class LooseHostnameVerifier implements HostnameVerifier {
	    public boolean verify(String hostname, SSLSession session) {
	        return true;
	    }
	}
}

```

```java:GeojsonProcesser.java
package citygml_lod1.components;

import java.awt.Shape;
import java.awt.geom.AffineTransform;
import java.awt.geom.Point2D;
import java.awt.geom.Rectangle2D;
import java.io.BufferedReader;
import java.io.File;
import java.io.FileReader;
import java.io.IOException;
import java.util.ArrayList;
import java.util.List;

import org.locationtech.jts.awt.ShapeReader;
import org.locationtech.jts.awt.ShapeWriter;
import org.locationtech.jts.geom.GeometryFactory;
import org.locationtech.jts.io.ParseException;
import org.locationtech.jts.io.geojson.GeoJsonReader;
import org.locationtech.jts.io.geojson.GeoJsonWriter;

import com.google.gson.JsonObject;
import com.mapbox.geojson.Feature;
import com.mapbox.geojson.FeatureCollection;
import com.mapbox.geojson.Geometry;
import com.mapbox.geojson.LineString;
import com.mapbox.geojson.MultiLineString;
import com.mapbox.geojson.MultiPoint;
import com.mapbox.geojson.MultiPolygon;
import com.mapbox.geojson.Point;
import com.mapbox.geojson.Polygon;

public class GeojsonProcesser {
	private List<Feature> list;
	private GeoJsonReader geomReader;
	private GeoJsonWriter geomWriter;
	private ShapeReader shapeReader;
	private ShapeWriter shapeWriter;
	private List<Shape> shapes;
	private Rectangle2D bounds;
	
	public GeojsonProcesser(FeatureCollection fc) {
		list=fc.features();
		geomReader=new GeoJsonReader();
		geomWriter=new GeoJsonWriter();
		shapeReader=new ShapeReader(new GeometryFactory());
		shapeWriter=new ShapeWriter();
		createShape();
	}
	
	public GeojsonProcesser(File f) throws IOException {
		this(FeatureCollection.fromJson(getString(f)));
	}
	
	public static String getString(File f) throws IOException{
		StringBuffer sb=new StringBuffer();
		String line=null;
		BufferedReader br=new BufferedReader(new FileReader(f));
		while((line=br.readLine())!=null){
			sb.append(line);
		}
		br.close();
		return sb.toString();
	}
	
	public FeatureCollection getFeatureCollection() {
		return FeatureCollection.fromFeatures(list);
	}
	
	public List<Feature> getFeatures() {
		return list;
	}
	
	public Feature getFeature(int i) {
		return list.get(i);
	}
	
	public JsonObject getProperty(int i){
		return list.get(i).properties();
	}
	
	public Geometry getGeometry(int i) {
		return list.get(i).geometry();
	}
	
	public Feature getFeature(Point2D p) {
		if(!bounds.contains(p))return null;
		int n=shapes.size();
		for(int i=0;i<n;i++) {
			Shape s=shapes.get(i);
			if(s.contains(p))return list.get(i);
		}
		return null;
	}
	
	public Feature getFeature(double x,double y) {
		if(!bounds.contains(x,y))return null;
		int n=shapes.size();
		for(int i=0;i<n;i++) {
			Shape s=shapes.get(i);
			if(s.contains(x,y))return list.get(i);
		}
		return null;
	}

	public List<Shape> getShapes(){
		return shapes;
	}
	
	public Shape getShape(int i) {
		return shapes.get(i);
	}
	
	public void updateShape() {
		int n=shapes.size();
		AffineTransform af=AffineTransform.getScaleInstance(1.0, 1.0);
		for(int i=0;i<n;i++) {
			Shape s=shapes.get(i);
			org.locationtech.jts.geom.Geometry g=shapeReader.read(s.getPathIterator(af));
			Geometry geo=reverseGeom(g);
			Feature f=list.get(i);
			Feature f2=Feature.fromGeometry(geo, f.properties());
			list.remove(i);
			list.add(i, f2);
		}
	}
	
	private void createShape(){
		bounds=null;
		shapes=new ArrayList<>();
		for(Feature f : list) {
			org.locationtech.jts.geom.Geometry g=transGeom(f.geometry());
			Shape sp=shapeWriter.toShape(g);
			shapes.add(sp);
			if(bounds==null) {
				bounds=sp.getBounds2D();
			}else {
				bounds.add(sp.getBounds2D());
			}
		}
	}
	
	private org.locationtech.jts.geom.Geometry transGeom(com.mapbox.geojson.Geometry geo){
		try {
			org.locationtech.jts.geom.Geometry geom=geomReader.read(geo.toJson());
			return geom;
		} catch (ParseException e) {
			e.printStackTrace();
			return null;
		}
	}
	
	private com.mapbox.geojson.Geometry reverseGeom(org.locationtech.jts.geom.Geometry geo){
		String json=geomWriter.write(geo);
		String type=geo.getGeometryType();
		if(type.equals(org.locationtech.jts.geom.Geometry.TYPENAME_LINESTRING)) {
			return LineString.fromJson(json);
		}else if(type.equals(org.locationtech.jts.geom.Geometry.TYPENAME_POINT)) {
			return Point.fromJson(json);
		}else if(type.equals(org.locationtech.jts.geom.Geometry.TYPENAME_POLYGON)) {
			return Polygon.fromJson(json);
		}else if(type.equals(org.locationtech.jts.geom.Geometry.TYPENAME_MULTILINESTRING)) {
			return MultiLineString.fromJson(json);
		}else if(type.equals(org.locationtech.jts.geom.Geometry.TYPENAME_MULTIPOINT)) {
			return MultiPoint.fromJson(json);
		}else if(type.equals(org.locationtech.jts.geom.Geometry.TYPENAME_MULTIPOLYGON)) {
			return MultiPolygon.fromJson(json);
		}
		return null;
	}
}
```

```java:TileApp.java

package citygml_lod1.app;

import java.awt.BorderLayout;
import java.awt.Font;
import java.awt.Shape;
import java.awt.datatransfer.DataFlavor;
import java.awt.datatransfer.Transferable;
import java.awt.datatransfer.UnsupportedFlavorException;
import java.awt.dnd.DnDConstants;
import java.awt.dnd.DropTarget;
import java.awt.dnd.DropTargetAdapter;
import java.awt.dnd.DropTargetDragEvent;
import java.awt.dnd.DropTargetDropEvent;
import java.awt.dnd.DropTargetListener;
import java.awt.event.WindowAdapter;
import java.awt.event.WindowEvent;
import java.awt.geom.AffineTransform;
import java.awt.geom.Rectangle2D;
import java.io.BufferedWriter;
import java.io.File;
import java.io.FileWriter;
import java.io.IOException;
import java.util.List;

import javax.imageio.ImageIO;
import javax.swing.JFrame;
import javax.swing.JLabel;
import javax.swing.JOptionPane;
import javax.swing.SwingUtilities;
import javax.swing.UIManager;

import citygml_lod1.components.GeojsonProcesser;
import citygml_lod1.components.WebTile;

public class TiletApp {
	private JFrame frame;
	private JLabel label;
	private String dsm="https://gio.pref.hyogo.lg.jp/tile/dsm/{z}/{y}/{x}.png";
	private String dem="https://gio.pref.hyogo.lg.jp/tile/dem/{z}/{y}/{x}.png";
	
	public TiletApp() {
		frame=new JFrame();
		frame.setTitle("DMgetApp");
		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.setSize(400,400);
		frame.setResizable(false);
		label=new JLabel("File Drop");
		label.setFont(new Font(Font.SANS_SERIF,Font.PLAIN,24));
		label.setHorizontalAlignment(JLabel.CENTER);
		label.setVerticalAlignment(JLabel.CENTER);
		frame.getContentPane().setLayout(new BorderLayout());
		frame.getContentPane().add(label,BorderLayout.CENTER);
		DropTargetListener dtl = new DropTargetAdapter() {
			  @Override public void dragOver(DropTargetDragEvent dtde) {
			    if (dtde.isDataFlavorSupported(DataFlavor.javaFileListFlavor)) {
			      dtde.acceptDrag(DnDConstants.ACTION_COPY);
			      return;
			    }
			    dtde.rejectDrag();
			  }
			  @SuppressWarnings("rawtypes")
			  @Override
			  public void drop(DropTargetDropEvent dtde) {
			    try {
			      if (dtde.isDataFlavorSupported(DataFlavor.javaFileListFlavor)) {
			        dtde.acceptDrop(DnDConstants.ACTION_COPY);
			        Transferable transferable = dtde.getTransferable();
			        List list = (List) transferable.getTransferData(DataFlavor.javaFileListFlavor);
			        for (Object o: list) {
			          if (o instanceof File) {
			            File f=(File) o;
			            if(f.getName().toLowerCase().endsWith(".geojson")) {
			            	proc(f);
			            }
			          }
			        }
			        dtde.dropComplete(true);
			        return;
			      }
			    } catch (UnsupportedFlavorException | IOException ex) {
			      ex.printStackTrace();
			    }
			    dtde.rejectDrop();
			  }
			};
			new DropTarget(label, DnDConstants.ACTION_COPY, dtl, true);
	}
	
	private void proc(File f) {
		Runnable r=new Runnable() {
			public void run() {
				try {
					setMessage("領域確認...");
					GeojsonProcesser gp = new GeojsonProcesser(f);
					Rectangle2D rect=null;
					for(Shape s : gp.getShapes()) {
						if(rect==null) {
							rect=s.getBounds2D();
						}else {
							rect.add(s.getBounds2D());
						}
					}
					setMessage("DSM取得...");
					WebTile tile=new WebTile(dsm,17,1);
					tile.create(5, rect);
					ImageIO.write(tile.getImage(), "png", new File(f.getAbsolutePath().replace(".geojson", "_dsm.png")));
					outTransform(tile.getTransform(),new File(f.getAbsolutePath().replace(".geojson", "_dsm.pgw")));
					setMessage("DEM取得...");
					tile=new WebTile(dem,17,1);
					tile.create(5, rect);
					ImageIO.write(tile.getImage(), "png", new File(f.getAbsolutePath().replace(".geojson", "_dem.png")));
					outTransform(tile.getTransform(),new File(f.getAbsolutePath().replace(".geojson", "_dem.pgw")));
					setMessage("終了...");
				} catch (IOException e) {
					e.printStackTrace();
				}

			}
		};
		new Thread(r).start();
	}
	
	private void outTransform(AffineTransform af,File path) throws IOException {
		BufferedWriter bw=new BufferedWriter(new FileWriter((path)));
		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.flush();
		bw.close();
	}
	
	private void setMessage(String args) {
		SwingUtilities.invokeLater(new Runnable() {
			public void run() {
				label.setText(args);
			}
		});
	}
	
	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);
		}
	}
	
	public static void main(String[] args) {
		TiletApp app=new TiletApp();
		app.frame.setLocationRelativeTo(null);
		app.frame.setVisible(true);
	}
}

```
TileAppを起動すると、ウィンドウが表示されます。

![015.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/55735/d8193c04-02ab-50be-4f47-9e05494090f7.jpeg)


Geojsonファイルをウィンドウ内にDropすると処理が開始され、しばらく（結構かも）待つと兵庫県高精度3次元データのDSM画像、DEM画像とワールドファイルが出力されます。

![005.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/55735/70f418c9-0892-4e94-bf2d-d8b12d3badc4.jpeg)

#3.Geojsonへの高さ属性の設定
出力されたDSM、DEMを使い、建物に地盤高（DEM）と表層高（DSM）の属性を追加します。
処理方法は、GeojsonからFeatureを逐次抽出し、ポリゴン（geometry.coordinates）内に点を発生させ、その点の表層高、地盤高をDSM、DEMから取得し、属性（properties）を追加するといった内容になります。
画像座標と地理座標（平面直角第5系）の変換は、DEM等を取得した際に出力したワールドファイルを使用します。ワールドファイルはアフィン変換行列を記述したファイルですので、これからAffineTransformインスタンスを生成し、その逆変換で地理座標を画像座標に変換します。
実装したクラスは以下のとおりです。
・GeoJsonSetDCHM.java：建物GeojsonのポリゴンにDEM、DSM属性を追加する。

```java:GeojsonSetDCHM.java

package citygml_lod1.app;

import java.awt.Shape;
import java.awt.geom.AffineTransform;
import java.awt.geom.NoninvertibleTransformException;
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.ArrayList;
import java.util.List;

import javax.imageio.ImageIO;

import citygml_lod1.components.GeojsonProcesser;

public class GeoJsonSetDCHM {

	public static void main(String[] args) {
		try {
			GeojsonProcesser geo=new GeojsonProcesser(new File("data/sannomiya.geojson"));
			BufferedImage dsm=ImageIO.read(new File("data/sannomiya_dsm.png"));
			BufferedImage dem=ImageIO.read(new File("data/sannomiya_dem.png"));
			AffineTransform af=loadTransform(new File("data/sannomiya_dsm.pgw"));
			try {
				af=af.createInverse();
			} catch (NoninvertibleTransformException e) {
				e.printStackTrace();
			}
			int n=geo.getShapes().size();
			for(int i=0;i<n;i++) {
				if(i%10000==0)System.out.println(i);
				Shape s=geo.getShape(i);
				Point2D p=getPoint(s);
				double dsmv=getZ(dsm,af,p);
				double demv=getZ(dem,af,p);
				geo.getProperty(i).addProperty("DEM", demv);
				geo.getProperty(i).addProperty("DCHM", dsmv-demv);
			}
			BufferedWriter bw=new BufferedWriter(new FileWriter(new File("data/sannomiya_gml.geojson")));
			bw.write(geo.getFeatureCollection().toJson());
			bw.close();
		} catch (IOException e) {
			e.printStackTrace();
		}		
	}
	
	private static double getZ(BufferedImage img,AffineTransform af,Point2D p) {
		Point2D px=af.transform(p, new Point2D.Double());
		int xx=(int)Math.round(px.getX());
		int yy=(int)Math.round(px.getY());
		if(xx>=0&&xx<img.getWidth()&&yy>=0&&yy<img.getHeight()) {
			return getVal(img.getRGB(xx, yy));
		}else {
			return 0.0;
		}
	}
	
	public static double getVal(int color){
		color=(color << 8) >> 8;
		if(color==8388608||color==-8388608){
			return Double.NaN;
		}else if(color<8388608){
			return color * 0.01;
		}else{
			return (color-16777216)*0.01;
		}
	}
	
	private static Point2D getPoint(Shape s) {
		Rectangle2D r=s.getBounds2D();
		double xx=0.0,yy=0.0;
		boolean flg=true;
		while(flg) {
			xx=Math.random()*r.getWidth()+r.getX();
			yy=Math.random()*r.getHeight()+r.getY();
			if(s.contains(xx, yy))flg=false;
		}
		return new Point2D.Double(xx,yy);
	}
	
	public static AffineTransform loadTransform(File path)throws IOException{
		BufferedReader br=new BufferedReader(new FileReader(path));
		List<Double> dd=new ArrayList<Double>();
		String line=null;
		while((line=br.readLine())!=null){
			double d=Double.parseDouble(line);
			dd.add(d);
		}
		br.close();
		double[] p=new double[dd.size()];
		for(int i=0;i<p.length;i++){
			p[i]=dd.get(i);
		}
		return new AffineTransform(p);
	}
}

```

上記処理を実行すると、建物に地盤高（DEM）と建物高さ（DCHM）が追加されたgeojsonファイルが出力されます。

#4.CityGML（LOD1）の出力
建物の平面形状（geometry.coordinates）と地盤高（properties.dem）、建物高さ（properties.dchm）からLOD1のbldg：境界面、bldg：屋根面、bldg:壁面、bldg;接地面を形成し、CityGML（LOD1）を生成します。
CityGMLの処理には[citygml4j](https://github.com/citygml4j/citygml4j)を使用しています。
実装したクラスは以下のとおりです。
・BldgCreator.java：建物Geojsonから、CityGML(LOD1)を生成するクラス。

```java:BldgCreator.java

package citygml_lod1.app;

import java.awt.Shape;
import java.awt.geom.AffineTransform;
import java.awt.geom.PathIterator;
import java.awt.geom.Point2D;
import java.io.File;
import java.io.IOException;
import java.util.ArrayList;
import java.util.List;

import org.citygml4j.CityGMLContext;
import org.citygml4j.builder.jaxb.CityGMLBuilder;
import org.citygml4j.builder.jaxb.CityGMLBuilderException;
import org.citygml4j.factory.DimensionMismatchException;
import org.citygml4j.factory.GMLGeometryFactory;
import org.citygml4j.geometry.BoundingBox;
import org.citygml4j.geometry.Point;
import org.citygml4j.model.citygml.CityGMLClass;
import org.citygml4j.model.citygml.building.AbstractBoundarySurface;
import org.citygml4j.model.citygml.building.BoundarySurfaceProperty;
import org.citygml4j.model.citygml.building.Building;
import org.citygml4j.model.citygml.building.GroundSurface;
import org.citygml4j.model.citygml.building.RoofSurface;
import org.citygml4j.model.citygml.building.WallSurface;
import org.citygml4j.model.citygml.core.CityModel;
import org.citygml4j.model.citygml.core.CityObjectMember;
import org.citygml4j.model.citygml.generics.MeasureAttribute;
import org.citygml4j.model.citygml.generics.StringAttribute;
import org.citygml4j.model.gml.basicTypes.Measure;
import org.citygml4j.model.gml.feature.BoundingShape;
import org.citygml4j.model.gml.geometry.aggregates.MultiSurface;
import org.citygml4j.model.gml.geometry.aggregates.MultiSurfaceProperty;
import org.citygml4j.model.gml.geometry.complexes.CompositeSurface;
import org.citygml4j.model.gml.geometry.primitives.Envelope;
import org.citygml4j.model.gml.geometry.primitives.Polygon;
import org.citygml4j.model.gml.geometry.primitives.Solid;
import org.citygml4j.model.gml.geometry.primitives.SolidProperty;
import org.citygml4j.model.gml.geometry.primitives.SurfaceProperty;
import org.citygml4j.model.gml.measures.Length;
import org.citygml4j.model.module.citygml.CityGMLVersion;
import org.citygml4j.util.bbox.BoundingBoxOptions;
import org.citygml4j.xml.io.CityGMLOutputFactory;
import org.citygml4j.xml.io.writer.CityGMLWriteException;
import org.citygml4j.xml.io.writer.CityGMLWriter;
import org.locationtech.jts.io.ParseException;

import com.google.gson.JsonObject;

import citygml_lod1.components.GeojsonProcesser;
import citygml_lod1.components.LonLatXY;

public class BldgCreator {
	private GeojsonProcesser geo;
	private int coordSys;

	public BldgCreator(File geojson,int cs) throws IOException, ParseException, CityGMLBuilderException {
		geo=new GeojsonProcesser(geojson);
		this.coordSys=cs;
	}
	
	public void create(File outfile) throws DimensionMismatchException, CityGMLWriteException, CityGMLBuilderException {
		CityGMLContext ctx = CityGMLContext.getInstance();
		CityGMLBuilder builder = ctx.createCityGMLBuilder();		
		CityModel cityModel = new CityModel();
		Point low=null,up=null;
		int n=geo.getShapes().size();
		for(int i=0;i<n;i++) {
			Shape s=geo.getShape(i);
			JsonObject p=geo.getProperty(i);
			double dchm=p.get("DCHM").getAsDouble();
			if(dchm<2.0)continue;
			double area=p.get("AREA").getAsDouble();
			if(area==0.0)continue;
			Building bldg=createBuilding(s,p);
			BoundingShape bs=bldg.calcBoundedBy(BoundingBoxOptions.defaults());
			BoundingBox bb=bs.getEnvelope().toBoundingBox();
			if(low==null) {
				low=bb.getLowerCorner();
			}else {
				Point ll=bb.getLowerCorner();
				low.setX(Math.min(low.getX(), ll.getX()));
				low.setY(Math.min(low.getY(), ll.getY()));
				low.setZ(Math.min(low.getZ(), ll.getZ()));
			}
			if(up==null) {
				up=bb.getUpperCorner();
			}else {
				Point ll=bb.getLowerCorner();
				up.setX(Math.max(up.getX(), ll.getX()));
				up.setY(Math.max(up.getY(), ll.getY()));
				up.setZ(Math.max(up.getZ(), ll.getZ()));
			}
//			cityModel.setBoundedBy(bldg.calcBoundedBy(BoundingBoxOptions.defaults()));
			cityModel.addCityObjectMember(new CityObjectMember(bldg));
		}
		BoundingBox bb=new BoundingBox();
		bb.setLowerCorner(low);
		bb.setUpperCorner(up);
		Envelope env=new BoundingShape(bb).getEnvelope();
		env.setSrsName("http://www.opengis.net/def/crs/EPSG/0/6673");
		cityModel.setBoundedBy(new BoundingShape(env));
		CityGMLOutputFactory out = builder.createCityGMLOutputFactory(CityGMLVersion.DEFAULT);
		CityGMLWriter writer = out.createCityGMLWriter(outfile, "UTF-8");
		writer.setPrefixes(CityGMLVersion.DEFAULT);
		writer.setSchemaLocations(CityGMLVersion.DEFAULT);
		writer.setIndentString("  ");
		writer.write(cityModel);
		writer.close();	
	}
	
	private Building createBuilding(Shape geo,JsonObject prop) throws DimensionMismatchException {
		GMLGeometryFactory geom = new GMLGeometryFactory();
		Building building = new Building();
		List<Point2D> pt=getPoints(geo);
		if(pt.size()<3)return null;
		double dchm=prop.get("DCHM").getAsDouble();
		double dem=prop.get("DEM").getAsDouble();
		Polygon ground = createGroundRoof(geom,pt,dem);
		List<Polygon> wall=createWall(geom,pt,dem,dem+dchm);
		Polygon roof = createGroundRoof(geom,pt,dchm+dem);
		StringAttribute sa=new StringAttribute();
        sa.setName("建物ID");
        sa.setValue(prop.get("gml_id").toString());        
        building.addGenericAttribute(sa);
        MeasureAttribute ma=new MeasureAttribute();
        Measure mo=new Measure(prop.get("PERIMETER").getAsDouble());
        mo.setUom("m");
        ma.setName("周辺長");
        ma.setValue(mo);
        building.addGenericAttribute(ma);
        ma=new MeasureAttribute();
        mo=new Measure(prop.get("AREA").getAsDouble());
        mo.setUom("m2");
        ma.setName("面積");
        ma.setValue(mo);
        building.addGenericAttribute(ma);
        sa=new StringAttribute();
        sa.setName("建物タイプ");
        sa.setValue(prop.get("type").toString());        
        building.addGenericAttribute(sa);
        Length ll=new Length();
        ll.setUom("m");
        ll.setValue(dchm);
        building.setMeasuredHeight(ll);
        building.setLod0FootPrint(new MultiSurfaceProperty(new MultiSurface(ground)));
        building.setLod0RoofEdge(new MultiSurfaceProperty(new MultiSurface(roof)));
		List<SurfaceProperty> surfaceMember = new ArrayList<>();
		surfaceMember.add(new SurfaceProperty(ground));
		for(Polygon pl : wall) {
			surfaceMember.add(new SurfaceProperty(pl));
		}
		surfaceMember.add(new SurfaceProperty(roof));
		
		/*
		surfaceMember.add(new SurfaceProperty('#' + ground.getId()));
		for(Polygon pl : wall) {
			surfaceMember.add(new SurfaceProperty('#' + pl.getId()));
		}
		surfaceMember.add(new SurfaceProperty('#' + roof.getId()));
		*/
		CompositeSurface compositeSurface = new CompositeSurface();
		compositeSurface.setSurfaceMember(surfaceMember);
		Solid solid = new Solid();
		solid.setSrsDimension(3);
		solid.setSrsName("http://www.opengis.net/def/crs/EPSG/0/6673");
		solid.setExterior(new SurfaceProperty(compositeSurface));
		building.setLod1Solid(new SolidProperty(solid));
		return building;
	}
	
	public BoundarySurfaceProperty createBoundarySurface(CityGMLClass type, Polygon geometry) {
		AbstractBoundarySurface boundarySurface = null;
		switch (type) {
		case BUILDING_WALL_SURFACE:
			boundarySurface = new WallSurface();
			break;
		case BUILDING_ROOF_SURFACE:
			boundarySurface = new RoofSurface();
			break;
		case BUILDING_GROUND_SURFACE:
			boundarySurface = new GroundSurface();
			break;
		default:
			break;
		}
		if (boundarySurface != null) {
			boundarySurface.setLod2MultiSurface(new MultiSurfaceProperty(new MultiSurface(geometry)));
			return new BoundarySurfaceProperty(boundarySurface);
		}
		return null;
	}
	
	private List<Polygon> createWall(GMLGeometryFactory geom,List<Point2D> l,double z1,double z2) throws DimensionMismatchException {
		List<Polygon> ret=new ArrayList<>();
		int n=l.size();
		for(int i=1;i<n;i++) {
			Point2D p1=l.get(i-1);
			Point2D p2=l.get(i);
			double[] p=new double[3*4];
			p[0]=p1.getX();
			p[1]=p1.getY();
			p[2]=z1;
			p[3]=p2.getX();
			p[4]=p2.getY();
			p[5]=z1;
			p[6]=p2.getX();
			p[7]=p2.getY();
			p[8]=z2;
			p[9]=p1.getX();
			p[10]=p1.getY();
			p[11]=z2;
			ret.add(geom.createLinearPolygon(p,3));
		}
		return ret;
	}
	
	private Polygon createGroundRoof(GMLGeometryFactory geom,List<Point2D> l,double z) throws DimensionMismatchException {
		int n=l.size();
		double[] p=new double[(n-1)*3];
		int it=0;
		for(int i=0;i<n-1;i++) {
			Point2D pt=l.get(i);
			p[it++]=pt.getX();
			p[it++]=pt.getY();
			p[it++]=z;
		}
		return geom.createLinearPolygon(p,3);
	}
	
	private List<Point2D> counterList(List<Point2D> l){
		List<Point2D> ret=new ArrayList<>();
		while(l.size()>0) {
			ret.add(l.remove(l.size()-1));
		}
		return ret;
	}
	
	private boolean isAntieClockWise(List<Point2D> p) {
		Point2D p1=p.get(0);
		Point2D p2=p.get(1);
		Point2D p3=p.get(2);
		double[] v1=new double[] {p2.getX()-p1.getX(),p2.getY()-p1.getY(),0};
		double[] v2=new double[] {p3.getX()-p1.getX(),p3.getY()-p1.getY(),0};
		double[] pd=crossProduct(v1,v2);
		if(pd[2]>0) {
			return true;
		}else {
			return false;
		}
	}
	
	
	private static double[] crossProduct(double[] a, double[] b) {  
	    double[] entries = new double[] {
	          a[1] * b[2] - a[2] * b[1],
	          a[2] * b[0] - a[0] * b[2],
	          a[0] * b[1] - a[1] * b[0]};
	    return entries;
	}
	
	private List<Point2D> getPoints(Shape sp){
		List<Point2D> ret=new ArrayList<>();
		PathIterator pi=sp.getPathIterator(AffineTransform.getScaleInstance(1.0, 1.0));
		double[] p=new double[6];
		while(!pi.isDone()){
			switch(pi.currentSegment(p)){
				case PathIterator.SEG_MOVETO:
					ret.add(new Point2D.Double(p[0],p[1]));
					break;
				case PathIterator.SEG_LINETO:
					ret.add(new Point2D.Double(p[0],p[1]));
					break;
//				case PathIterator.SEG_QUADTO:
//					ret.add(new Point2D.Double(p[2],p[3]));
//					break;
//				case PathIterator.SEG_CUBICTO:
//					ret.add(new Point2D.Double(p[4],p[5]));
//					break;
				case PathIterator.SEG_CLOSE:
					Point2D c=ret.get(0);
					ret.add(new Point2D.Double(c.getX(),c.getY()));
					break;
			}
			pi.next();
		}
		if(isAntieClockWise(ret)) {
			return transLonlat(ret);
		}else {
			return transLonlat(counterList(ret));
		}
	}
	
	private List<Point2D> transLonlat(List<Point2D> l){
		List<Point2D> ret=new ArrayList<>();
		for(Point2D p : l) {
			Point2D pt=LonLatXY.xyToLonlat(coordSys, p.getX(), p.getY());
			ret.add(new Point2D.Double(pt.getY(),pt.getX()));
		}
		return ret;
	}
	
	public static void main(String[] args) {
		File f=new File("data/sannomiya_gml.geojson");
		File o=new File("data/sannomiya.gml");
		try {
			BldgCreator bc=new BldgCreator(f,5);
			bc.create(o);
		} catch (IOException | ParseException | CityGMLBuilderException | DimensionMismatchException | CityGMLWriteException e) {
			e.printStackTrace();
		}
	}
}

```

上記処理を実行すると、三宮周辺のCityGML（sannomiya.gml）が出力されます。

生成したCityGMLを[FZKViewer](https://www.iai.kit.edu/1302.php)で表示すると以下のようになりました。
検査した訳ではないですが、まあ三宮っぽくはあります。

![011.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/55735/6b384b78-8cb2-f8ee-e5f5-afa406c20e61.jpeg)
![012.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/55735/883354b3-bd4d-8716-f8e4-f5ff389fe6db.jpeg)

#5.まとめ
Javaで、基盤地図情報とDEM、DSMデータから、CityGML（LOD1）を生成してみました。
mapbox-gl等にのせる場合は、Mapbox Vector Tileに変換して使うと良い感じです。

https://qiita.com/t-mat/items/a57f44c46d9e7f0aac56

LOD1は、2D建物ポリゴンに建物高さ属性を追加したものと、ほとんど同じですので、DSMとDEMがあれば、比較的簡単に作れる感じがします。
課題は、基盤地図情報の建物ポリゴンと実際の建物状況（航空写真など）が異なる場合があることや、DSMデータの入手が難しいことが挙げられます。