Commit f0f212e1 authored by Matthias Betz's avatar Matthias Betz
Browse files

add wkt polygon fitting

parent fe063002
......@@ -12,6 +12,7 @@ import java.util.HashSet;
import java.util.List;
import java.util.Map;
import java.util.Set;
import org.citygml4j.core.model.core.CityModel;
import org.citygml4j.xml.CityGMLContextException;
import org.citygml4j.xml.reader.CityGMLReadException;
......@@ -21,10 +22,19 @@ import org.geotools.data.DataStoreFinder;
import org.geotools.data.FeatureSource;
import org.geotools.feature.FeatureCollection;
import org.geotools.feature.FeatureIterator;
import org.locationtech.jts.geom.CoordinateFilter;
import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.geom.MultiPolygon;
import org.locationtech.jts.geom.Polygon;
import org.locationtech.jts.io.ParseException;
import org.locationtech.jts.io.WKTReader;
import org.locationtech.proj4j.BasicCoordinateTransform;
import org.locationtech.proj4j.CRSFactory;
import org.locationtech.proj4j.CoordinateReferenceSystem;
import org.locationtech.proj4j.ProjCoordinate;
import org.opengis.feature.simple.SimpleFeature;
import org.opengis.feature.simple.SimpleFeatureType;
import de.hft.stuttgart.citygml.green.osm.GreenArea;
import de.hft.stuttgart.citygml.green.osm.GreenEnricher;
import de.hft.stuttgart.citygml.green.osm.LandUseArea;
......@@ -36,6 +46,12 @@ import jakarta.xml.bind.JAXBException;
public class AlkisGreenEnricher
{
private static final CRSFactory CRS_FACTORY = new CRSFactory();
private static final String WGS_84 = "EPSG:4326";
private static final CoordinateReferenceSystem wgs84 = CRS_FACTORY.createFromName(WGS_84);
private static final CoordinateReferenceSystem targetCrs = CRS_FACTORY.createFromName("EPSG:25832");
private static final BasicCoordinateTransform bct = new BasicCoordinateTransform(wgs84, targetCrs);
private static Set<String> greenAreaTypes = new HashSet<>();
private static Set<String> roadAreaTypes = new HashSet<>();
......@@ -55,45 +71,69 @@ public class AlkisGreenEnricher
public static void main(String[] args)
throws IOException, InterruptedException, CityGMLContextException, CityGMLReadException, JAXBException,
CityGMLWriteException {
CityGMLWriteException, ParseException {
System.out.println("Reading CityGML file");
Path inFile = Paths.get(args[0]);
Path alkisDataPath = Paths.get(args[1]); // shp, shx and dbf should be present.
Path baumKatasterPath = Paths.get(args[2]);
String outputSuffix = args[3];
Polygon wktPolygon = null;
if (args.length == 5) {
WKTReader wktReader = new WKTReader();
Geometry wktGeometry = wktReader.read(args[4]);
if (!(wktGeometry instanceof Polygon)) {
throw new IllegalArgumentException("WKT is not a polygon");
}
wktPolygon = (Polygon) wktGeometry;
wktPolygon.apply((CoordinateFilter) c -> {
ProjCoordinate p1 = new ProjCoordinate();
p1.x = c.x;
p1.y = c.y;
ProjCoordinate p2 = new ProjCoordinate();
bct.transform(p1, p2);
c.x = p2.x;
c.y = p2.y;
});
wktPolygon.geometryChanged();
}
CityModel cityModel = GreenEnricher.readCityGml(inFile);
GreenEnricher.createTransformers(cityModel);
OsmData osmData = new OsmData();
String boundingBoxString = GreenEnricher.extractAndConvertBoundingBox(cityModel, osmData);
String boundingBoxBasename = boundingBoxString.replace(",", "__").replace('.', '_');
Path osmCache = Paths.get("data", "cache", "osm_response_" + boundingBoxBasename + ".xml");
if (!Files.exists(osmCache)) {
System.out.println("Downloading OSM data for " + boundingBoxString);
HttpResponse<String> response = GreenEnricher.getOsmData(boundingBoxString);
Files.write(osmCache, response.body().getBytes(StandardCharsets.UTF_8));
}
String osmResponse = Files.readString(osmCache);
System.out.println("Parsing OSM response");
GreenEnricher.parseOsmResponse(osmResponse, osmData);
// String boundingBoxString = GreenEnricher.extractAndConvertBoundingBox(cityModel, osmData);
// String boundingBoxBasename = boundingBoxString.replace(",", "__").replace('.', '_');
// Path osmCache = Paths.get("data", "cache", "osm_response_" + boundingBoxBasename + ".xml");
// if (!Files.exists(osmCache)) {
// System.out.println("Downloading OSM data for " + boundingBoxString);
// HttpResponse<String> response = GreenEnricher.getOsmData(boundingBoxString);
// Files.write(osmCache, response.body().getBytes(StandardCharsets.UTF_8));
// }
// String osmResponse = Files.readString(osmCache);
// System.out.println("Parsing OSM response");
// GreenEnricher.parseOsmResponse(osmResponse, osmData);
// ignore green areas from osm
osmData.getGreenAreas().clear();
// osmData.getGreenAreas().clear();
parseAlkisData(osmData, alkisDataPath);
System.out.println("Fit data in bounding box");
GreenEnricher.fitToBoundingBox(osmData);
if (wktPolygon != null) {
GreenEnricher.fitToPolygon(osmData, wktPolygon);
} else {
GreenEnricher.fitToBoundingBox(osmData);
}
GreenEnricher.convertGreenAreasToCityGML(cityModel, osmData.getGreenAreas());
GreenEnricher.convertWaterAreasToCityGML(cityModel, osmData);
GreenEnricher.convertRoadAreasToCityGML(cityModel, osmData);
GreenEnricher.converLandUseAreasToCityGML(cityModel, osmData);
GreenEnricher.convertLandUseAreasToCityGML(cityModel, osmData);
TreeUtils.insertTrees(cityModel, osmData, baumKatasterPath);
TreeUtils.insertTrees(cityModel, osmData, baumKatasterPath, wktPolygon);
GreenEnricher.clampToGround(cityModel);
......
......@@ -18,7 +18,7 @@ import java.util.Set;
import java.util.UUID;
import java.util.regex.Matcher;
import java.util.regex.Pattern;
import javax.xml.parsers.ParserConfigurationException;
import org.citygml4j.core.model.CityGMLVersion;
import org.citygml4j.core.model.building.Building;
import org.citygml4j.core.model.core.AbstractCityObject;
......@@ -40,12 +40,15 @@ import org.citygml4j.xml.writer.CityGMLOutputFactory;
import org.citygml4j.xml.writer.CityGMLWriteException;
import org.citygml4j.xml.writer.CityGMLWriter;
import org.locationtech.jts.geom.Coordinate;
import org.locationtech.jts.geom.CoordinateFilter;
import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.geom.GeometryFactory;
import org.locationtech.jts.geom.LineString;
import org.locationtech.jts.geom.MultiPolygon;
import org.locationtech.jts.geom.Point;
import org.locationtech.jts.geom.Polygon;
import org.locationtech.jts.io.ParseException;
import org.locationtech.jts.io.WKTReader;
import org.locationtech.proj4j.BasicCoordinateTransform;
import org.locationtech.proj4j.CRSFactory;
import org.locationtech.proj4j.CoordinateReferenceSystem;
......@@ -60,6 +63,7 @@ import org.xmlobjects.gml.model.geometry.aggregates.MultiSurfaceProperty;
import org.xmlobjects.gml.model.geometry.primitives.AbstractRingProperty;
import org.xmlobjects.gml.model.geometry.primitives.LinearRing;
import org.xmlobjects.gml.model.geometry.primitives.SurfaceProperty;
import de.hft.stuttgart.citygml.green.osm.jaxb.OSM;
import de.hft.stuttgart.citygml.green.osm.jaxb.OsmMember;
import de.hft.stuttgart.citygml.green.osm.jaxb.OsmNode;
......@@ -73,7 +77,6 @@ import jakarta.xml.bind.JAXBException;
public class GreenEnricher
{
private static final int BOUNDING_BOX_INCREASE_IN_M = 100;
// in degrees
......@@ -107,17 +110,38 @@ public class GreenEnricher
private static final URI OVERPASS_API_URI = URI.create("http://www.overpass-api.de/api/interpreter");
private static final CRSFactory CRS_FACTORY = new CRSFactory();
private static CoordinateReferenceSystem sourceCRS;
private static CoordinateReferenceSystem targetCRS = CRS_FACTORY.createFromName("EPSG:4326");
private static CoordinateReferenceSystem wgs84 = CRS_FACTORY.createFromName("EPSG:4326");
private static BasicCoordinateTransform transform;
private static BasicCoordinateTransform backTransform;
public static final GeometryFactory GEOM_FACTORY = new GeometryFactory();
public static void main(String[] args) throws IOException, CityGMLContextException, CityGMLReadException,
InterruptedException, ParserConfigurationException, CityGMLWriteException, JAXBException {
InterruptedException, CityGMLWriteException, JAXBException, ParseException {
System.out.println("Reading CityGML file");
Path inFile = Paths.get(args[0]);
Path baumKatasterPath = Paths.get(args[1]);
String outputSuffix = args[2];
Polygon wktPolygon = null;
if (args.length == 4) {
WKTReader wktReader = new WKTReader();
Geometry wktGeometry = wktReader.read(args[3]);
if (!(wktGeometry instanceof Polygon)) {
throw new IllegalArgumentException("WKT is not a polygon");
}
wktPolygon = (Polygon) wktGeometry;
CoordinateReferenceSystem utm32 = CRS_FACTORY.createFromName("EPSG:25832");
BasicCoordinateTransform bct = new BasicCoordinateTransform(wgs84, utm32);
wktPolygon.apply((CoordinateFilter) c -> {
ProjCoordinate p1 = new ProjCoordinate();
p1.x = c.x;
p1.y = c.y;
ProjCoordinate p2 = new ProjCoordinate();
bct.transform(p1, p2);
c.x = p2.x;
c.y = p2.y;
});
wktPolygon.geometryChanged();
}
CityModel cityModel = readCityGml(inFile);
createTransformers(cityModel);
......@@ -137,7 +161,11 @@ public class GreenEnricher
parseOsmResponse(osmResponse, osmData);
System.out.println("Fit data in bounding box");
fitToBoundingBox(osmData);
if (wktPolygon != null) {
fitToPolygon(osmData, wktPolygon);
} else {
fitToBoundingBox(osmData);
}
System.out.println("Filter intersecting areas");
List<GreenArea> greenAreas = osmData.getGreenAreas();
......@@ -147,7 +175,7 @@ public class GreenEnricher
convertWaterAreasToCityGML(cityModel, osmData);
// trees
TreeUtils.insertTrees(cityModel, osmData, baumKatasterPath);
TreeUtils.insertTrees(cityModel, osmData, baumKatasterPath, wktPolygon);
clampToGround(cityModel);
......@@ -161,11 +189,27 @@ public class GreenEnricher
}
public static void fitToPolygon(OsmData osmData, Polygon wktPolygon) {
fitGreenAreasToBounds(wktPolygon, osmData.getGreenAreas());
fitRoadAreasToBounds(wktPolygon, osmData.getRoadAreas());
fitLandUseAreasToBounds(wktPolygon, osmData.getLandUseAreas());
fitWaterAreasToBounds(wktPolygon, osmData.getWaterAreas());
for (Iterator<TreePoint> iterator = osmData.getTreePoints().iterator(); iterator.hasNext();) {
TreePoint tp = iterator.next();
if (!osmData.getBoundingBox().contains(tp.getPoint())) {
iterator.remove();
}
}
}
public static void createTransformers(CityModel cityModel) {
String epsgCode = extractEpsgCode(cityModel);
sourceCRS = CRS_FACTORY.createFromName(epsgCode);
transform = new BasicCoordinateTransform(sourceCRS, targetCRS);
backTransform = new BasicCoordinateTransform(targetCRS, sourceCRS);
transform = new BasicCoordinateTransform(sourceCRS, wgs84);
backTransform = new BasicCoordinateTransform(wgs84, sourceCRS);
}
......@@ -202,11 +246,10 @@ public class GreenEnricher
}
public static void fitToBoundingBox(OsmData osmData) {
fitGreenAreas(osmData, osmData.getGreenAreas());
fitRoadAreas(osmData, osmData.getRoadAreas());
fitLandUseAreas(osmData, osmData.getLandUseAreas());
clipWaterAreasToBoundingBox(osmData);
fitGreenAreasToBounds(osmData.getBoundingBox(), osmData.getGreenAreas());
fitRoadAreasToBounds(osmData.getBoundingBox(), osmData.getRoadAreas());
fitLandUseAreasToBounds(osmData.getBoundingBox(), osmData.getLandUseAreas());
fitWaterAreasToBounds(osmData.getBoundingBox(), osmData.getWaterAreas());
for (Iterator<TreePoint> iterator = osmData.getTreePoints().iterator(); iterator.hasNext();) {
TreePoint tp = iterator.next();
......@@ -217,11 +260,11 @@ public class GreenEnricher
}
private static void fitLandUseAreas(OsmData osmData, List<LandUseArea> landUseAreas) {
private static void fitLandUseAreasToBounds(Geometry bounds, List<LandUseArea> landUseAreas) {
List<LandUseArea> newLandUseAreas = new ArrayList<>();
for (LandUseArea landUseArea : landUseAreas) {
Polygon area = landUseArea.getArea();
Geometry intersection = area.intersection(osmData.getBoundingBox());
Geometry intersection = area.intersection(bounds);
if (intersection instanceof MultiPolygon multi) {
Polygon poly1 = (Polygon) multi.getGeometryN(0);
landUseArea.setArea(poly1);
......@@ -238,11 +281,11 @@ public class GreenEnricher
}
private static void fitRoadAreas(OsmData osmData, List<RoadArea> roadAreas) {
private static void fitRoadAreasToBounds(Geometry bounds, List<RoadArea> roadAreas) {
List<RoadArea> newRoadAreas = new ArrayList<>();
for (RoadArea roadArea : roadAreas) {
Polygon area = roadArea.getArea();
Geometry intersection = area.intersection(osmData.getBoundingBox());
Geometry intersection = area.intersection(bounds);
if (intersection instanceof MultiPolygon multi) {
Polygon poly1 = (Polygon) multi.getGeometryN(0);
roadArea.setArea(poly1);
......@@ -259,11 +302,11 @@ public class GreenEnricher
}
private static void fitGreenAreas(OsmData osmData, List<GreenArea> greenAreas) {
private static void fitGreenAreasToBounds(Geometry bounds, List<GreenArea> greenAreas) {
List<GreenArea> newGreenAreas = new ArrayList<>();
for (GreenArea greenArea : greenAreas) {
Polygon area = greenArea.getArea();
Geometry intersection = area.intersection(osmData.getBoundingBox());
Geometry intersection = area.intersection(bounds);
if (intersection instanceof MultiPolygon multi) {
Polygon poly1 = (Polygon) multi.getGeometryN(0);
greenArea.setArea(poly1);
......@@ -279,11 +322,11 @@ public class GreenEnricher
greenAreas.addAll(newGreenAreas);
}
private static void clipWaterAreasToBoundingBox(OsmData osmData) {
private static void fitWaterAreasToBounds(Geometry bounds, List<WaterArea> waterAreas) {
List<WaterArea> newWaterAreas = new ArrayList<>();
for (WaterArea waterArea : osmData.getWaterAreas()) {
for (WaterArea waterArea : waterAreas) {
Polygon area = waterArea.getArea();
Geometry intersection = area.intersection(osmData.getBoundingBox());
Geometry intersection = area.intersection(bounds);
if (intersection instanceof MultiPolygon multi) {
Polygon poly1 = (Polygon) multi.getGeometryN(0);
waterArea.setArea(poly1);
......@@ -296,6 +339,7 @@ public class GreenEnricher
waterArea.setArea((Polygon) intersection);
}
}
waterAreas.addAll(newWaterAreas);
}
private static void removeDuplicateAreas(List<GreenArea> greenAreas) {
......@@ -808,7 +852,7 @@ public class GreenEnricher
}
public static void converLandUseAreasToCityGML(CityModel cityModel, OsmData osmData) {
public static void convertLandUseAreasToCityGML(CityModel cityModel, OsmData osmData) {
for (LandUseArea landUseArea : osmData.getLandUseAreas()) {
LandUse landUse = new LandUse();
org.xmlobjects.gml.model.geometry.primitives.Polygon poly = convertToCityGmlPoly(landUseArea.getArea());
......
......@@ -8,26 +8,31 @@ import org.citygml4j.core.model.core.AbstractCityObjectProperty;
import org.citygml4j.core.model.core.CityModel;
import org.citygml4j.core.model.vegetation.SolitaryVegetationObject;
import org.locationtech.jts.geom.Coordinate;
import org.locationtech.jts.geom.GeometryFactory;
import org.locationtech.jts.geom.Polygon;
import org.xmlobjects.gml.model.basictypes.Code;
import org.xmlobjects.gml.model.geometry.aggregates.MultiSurface;
import org.xmlobjects.gml.model.geometry.aggregates.MultiSurfaceProperty;
public class TreeUtils {
public class TreeUtils
{
private static final GeometryFactory FACTORY = new GeometryFactory();
public static void insertTrees(CityModel cityModel, OsmData osmData, Path baumKatasterPath) throws IOException {
public static void insertTrees(CityModel cityModel, OsmData osmData, Path baumKatasterPath, Polygon wktPolygon)
throws IOException {
TreeKatasterData katasterData = TreeKatasterData.parseTreeKatasterData(baumKatasterPath);
generateTreesFromKataster(cityModel, katasterData);
generateTreesFromKataster(cityModel, katasterData, osmData, wktPolygon);
// TreeKatasterData katasterData2 = TreeKatasterData.parseTreeKatasterData(Paths.get("data", "Trees_realistic_20240201.shp"));
// generateTreesFromKataster(cityModel, katasterData2);
// TreeKatasterData katasterData2 =
// TreeKatasterData.parseTreeKatasterData(Paths.get("data",
// "Trees_realistic_20240201.shp"));
// generateTreesFromKataster(cityModel, katasterData2);
// All kataster trees are taken, osm trees are removed
filterDuplicateTreesFromOSM(osmData, katasterData);
// filterDuplicateTreesFromOSM(osmData, katasterData2);
// filterDuplicateTreesFromOSM(osmData, katasterData);
// filterDuplicateTreesFromOSM(osmData, katasterData2);
generateTreesFromOSM(cityModel, osmData);
// generateTreesFromOSM(cityModel, osmData, wktPolygon);
}
......@@ -43,14 +48,19 @@ public class TreeUtils
}
}
private static void generateTreesFromOSM(CityModel cityModel, OsmData osmData) {
private static void generateTreesFromOSM(CityModel cityModel, OsmData osmData, Polygon wktPolygon) {
for (TreePoint tp : osmData.getTreePoints()) {
Coordinate coordinate = tp.getPoint().getCoordinate();
if ((wktPolygon != null && !wktPolygon.contains(FACTORY.createPoint(coordinate)))
|| !osmData.getBoundingBox().contains(FACTORY.createPoint(coordinate))) {
// ignore trees outside bounds
continue;
}
// standard tree
double trunkRadius = 0.89 / (2 * Math.PI);
double trunkHeight = TreeKatasterData.TRUNK_PERCENTAGE * 11.46;
double crownHeight = TreeKatasterData.CROWN_PERCENTAGE * 11.46;
double crownRadius = 8 / 2d;
Coordinate coordinate = tp.getPoint().getCoordinate();
if (Double.isNaN(coordinate.z)) {
coordinate.z = 0;
}
......@@ -64,14 +74,21 @@ public class TreeUtils
}
}
private static void generateTreesFromKataster(CityModel cityModel, TreeKatasterData katasterData) {
private static void generateTreesFromKataster(CityModel cityModel, TreeKatasterData katasterData, OsmData osmData,
Polygon wktPolygon) {
for (Tree tree : katasterData.getTrees()) {
// check if tree is in area
Coordinate coordinate = tree.getPoint().getCoordinate();
if ((wktPolygon != null && !wktPolygon.contains(FACTORY.createPoint(coordinate)))
|| !osmData.getBoundingBox().contains(FACTORY.createPoint(coordinate))) {
continue;
}
double trunkRadius = tree.getTrunkRadius();
double trunkHeight = tree.getTrunkHeight();
double crownHeight = tree.getCrownHeight();
double crownRadius = tree.getCrownRadius();
Coordinate coordinate = tree.getPoint().getCoordinate();
if (Double.isNaN(coordinate.z)) {
coordinate.z = 0;
}
......
......@@ -21,7 +21,8 @@ class AlkisGreenEnricherTest
"data/tn_09663/Nutzung.shp", // ALKIS Data
//NOTE: From https://transfer.hft-stuttgart.de/gitlab/circulargreensimcity/circulargreensimcity/-/issues/26#note_6544
"data/Trees/Trees_realisticScenario_20240201.shp", // Added trees, in Baumkatasterformat,
"alkis_test" // Output GML suffix
"alkis_test", // Output GML suffix
"POLYGON((9.946874 49.803172, 9.949785 49.803205, 9.952122 49.803312, 9.953435 49.803089, 9.956488 49.803150, 9.955465 49.802374, 9.953216 49.800916, 9.952371 49.800676, 9.951810 49.800509, 9.949172 49.800829, 9.946792 49.800796, 9.946874 49.803172))"
};
AlkisGreenEnricher.main(args);
assertTrue(Files.exists(outputGML));
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment