From b96d17f3eb275a97109cc160db1cfd3731a01d31 Mon Sep 17 00:00:00 2001 From: Mathieu Baudier Date: Sat, 30 Sep 2023 11:01:12 +0200 Subject: [PATCH] Work on GPX polygons clean up --- .../app/geo/{GeoJSon.java => GeoJson.java} | 93 +++++++++------ .../src/org/argeo/app/geo/GeoUtils.java | 111 +++++++++++++++++- .../src/org/argeo/app/geo/GpxUtils.java | 14 ++- .../argeo/app/geo/http/WfsHttpHandler.java | 10 +- 4 files changed, 177 insertions(+), 51 deletions(-) rename org.argeo.app.geo/src/org/argeo/app/geo/{GeoJSon.java => GeoJson.java} (60%) diff --git a/org.argeo.app.geo/src/org/argeo/app/geo/GeoJSon.java b/org.argeo.app.geo/src/org/argeo/app/geo/GeoJson.java similarity index 60% rename from org.argeo.app.geo/src/org/argeo/app/geo/GeoJSon.java rename to org.argeo.app.geo/src/org/argeo/app/geo/GeoJson.java index 5d5b99f..c1ba1d3 100644 --- a/org.argeo.app.geo/src/org/argeo/app/geo/GeoJSon.java +++ b/org.argeo.app.geo/src/org/argeo/app/geo/GeoJson.java @@ -3,8 +3,10 @@ package org.argeo.app.geo; import org.locationtech.jts.geom.Coordinate; import org.locationtech.jts.geom.Envelope; import org.locationtech.jts.geom.Geometry; +import org.locationtech.jts.geom.GeometryCollection; import org.locationtech.jts.geom.LineString; import org.locationtech.jts.geom.LinearRing; +import org.locationtech.jts.geom.MultiPoint; import org.locationtech.jts.geom.Point; import org.locationtech.jts.geom.Polygon; @@ -17,13 +19,16 @@ import jakarta.json.stream.JsonGenerator; * * @see https://datatracker.ietf.org/doc/html/rfc7946 */ -public class GeoJSon { +public class GeoJson { public final static String POINT_TYPE = "Point"; + public final static String MULTI_POINT_TYPE = "MultiPoint"; public final static String LINE_STRING_TYPE = "LineString"; public final static String POLYGON_TYPE = "Polygon"; + public final static String GEOMETRY_COLLECTION_TYPE = "GeometryCollection"; public final static String TYPE = "type"; public final static String GEOMETRY = "geometry"; + public final static String GEOMETRIES = "geometries"; public final static String COORDINATES = "coordinates"; public final static String BBOX = "bbox"; public final static String PROPERTIES = "properties"; @@ -32,65 +37,81 @@ public class GeoJSon { * WRITE */ /** Writes a {@link Geometry} as GeoJSON. */ - public static void writeGeometry(JsonGenerator generator, Geometry geometry) { + public static void writeGeometry(JsonGenerator g, Geometry geometry) { if (geometry instanceof Point point) { - generator.write(TYPE, POINT_TYPE); - generator.writeStartArray(COORDINATES); - writeCoordinate(generator, point.getCoordinate()); - generator.writeEnd();// coordinates array + g.write(TYPE, POINT_TYPE); + g.writeStartArray(COORDINATES); + writeCoordinate(g, point.getCoordinate()); + g.writeEnd();// coordinates array + } else if (geometry instanceof MultiPoint multiPoint) { + g.write(TYPE, MULTI_POINT_TYPE); + g.writeStartArray(COORDINATES); + writeCoordinates(g, multiPoint.getCoordinates()); + g.writeEnd();// coordinates array } else if (geometry instanceof LineString lineString) { - generator.write(TYPE, LINE_STRING_TYPE); - generator.writeStartArray(COORDINATES); - writeCoordinates(generator, lineString.getCoordinates()); - generator.writeEnd();// coordinates array + g.write(TYPE, LINE_STRING_TYPE); + g.writeStartArray(COORDINATES); + writeCoordinates(g, lineString.getCoordinates()); + g.writeEnd();// coordinates array } else if (geometry instanceof Polygon polygon) { - generator.write(TYPE, POLYGON_TYPE); - generator.writeStartArray(COORDINATES); + g.write(TYPE, POLYGON_TYPE); + g.writeStartArray(COORDINATES); LinearRing exteriorRing = polygon.getExteriorRing(); - generator.writeStartArray(); - writeCoordinates(generator, exteriorRing.getCoordinates()); - generator.writeEnd(); + g.writeStartArray(); + writeCoordinates(g, exteriorRing.getCoordinates()); + g.writeEnd(); for (int i = 0; i < polygon.getNumInteriorRing(); i++) { LinearRing interiorRing = polygon.getInteriorRingN(i); // TODO verify that holes are clockwise - generator.writeStartArray(); - writeCoordinates(generator, interiorRing.getCoordinates()); - generator.writeEnd(); + g.writeStartArray(); + writeCoordinates(g, interiorRing.getCoordinates()); + g.writeEnd(); } - generator.writeEnd();// coordinates array + g.writeEnd();// coordinates array + } else if (geometry instanceof GeometryCollection geometryCollection) {// must be last + g.write(TYPE, GEOMETRY_COLLECTION_TYPE); + g.writeStartArray(GEOMETRIES); + for (int i = 0; i < geometryCollection.getNumGeometries(); i++) { + g.writeStartObject(); + writeGeometry(g, geometryCollection.getGeometryN(i)); + g.writeEnd();// geometry object + } + g.writeEnd();// geometries array + } else { + throw new IllegalArgumentException(geometry.getClass() + " is not supported."); } } /** Writes a sequence of coordinates [[lat,lon],[lat,lon]] */ - public static void writeCoordinates(JsonGenerator generator, Coordinate[] coordinates) { + public static void writeCoordinates(JsonGenerator g, Coordinate[] coordinates) { for (Coordinate coordinate : coordinates) { - generator.writeStartArray(); - writeCoordinate(generator, coordinate); - generator.writeEnd(); + g.writeStartArray(); + writeCoordinate(g, coordinate); + g.writeEnd(); } } /** Writes a pair of coordinates [lat,lon]. */ - public static void writeCoordinate(JsonGenerator generator, Coordinate coordinate) { - generator.write(coordinate.getX()); - generator.write(coordinate.getY()); + public static void writeCoordinate(JsonGenerator g, Coordinate coordinate) { + g.write(coordinate.getX()); + g.write(coordinate.getY()); double z = coordinate.getZ(); if (!Double.isNaN(z)) { - generator.write(z); + g.write(z); } } /** * Writes the {@link Envelope} of a {@link Geometry} as a bbox GeoJSON object. */ - public static void writeBBox(JsonGenerator generator, Geometry geometry) { - generator.writeStartArray(BBOX); + public static void writeBBox(JsonGenerator g, Geometry geometry) { + g.writeStartArray(BBOX); Envelope envelope = geometry.getEnvelopeInternal(); - generator.write(envelope.getMinX()); - generator.write(envelope.getMinY()); - generator.write(envelope.getMaxX()); - generator.write(envelope.getMaxY()); - generator.writeEnd(); + g.write(envelope.getMinX()); + g.write(envelope.getMinY()); + g.write(envelope.getMaxX()); + g.write(envelope.getMaxY()); + g.writeEnd(); } /* @@ -125,7 +146,7 @@ public class GeoJSon { throw new IllegalArgumentException("Unexpected value: " + type); }; // res.normalize(); - return (T)res; + return (T) res; } /** Reads a coordinate sequence [[lat,lon],[lat,lon]]. */ @@ -143,7 +164,7 @@ public class GeoJSon { } /** singleton */ - private GeoJSon() { + private GeoJson() { } } diff --git a/org.argeo.app.geo/src/org/argeo/app/geo/GeoUtils.java b/org.argeo.app.geo/src/org/argeo/app/geo/GeoUtils.java index 8f97129..2c0bb39 100644 --- a/org.argeo.app.geo/src/org/argeo/app/geo/GeoUtils.java +++ b/org.argeo.app.geo/src/org/argeo/app/geo/GeoUtils.java @@ -8,6 +8,7 @@ import java.util.ArrayList; import java.util.HashMap; import java.util.List; import java.util.Map; +import java.util.Objects; import javax.measure.Quantity; import javax.measure.quantity.Area; @@ -44,9 +45,13 @@ import org.geotools.styling.css.CssParser; import org.geotools.styling.css.CssTranslator; import org.geotools.styling.css.Stylesheet; import org.geotools.xml.styling.SLDTransformer; +import org.locationtech.jts.algorithm.hull.ConcaveHull; import org.locationtech.jts.geom.Coordinate; +import org.locationtech.jts.geom.Geometry; import org.locationtech.jts.geom.Point; import org.locationtech.jts.geom.Polygon; +import org.locationtech.jts.geom.util.GeometryFixer; +import org.locationtech.jts.operation.polygonize.Polygonizer; import org.opengis.feature.simple.SimpleFeature; import org.opengis.feature.simple.SimpleFeatureType; import org.opengis.filter.Filter; @@ -68,17 +73,45 @@ public class GeoUtils { /** In square meters. */ public static Quantity calcArea(SimpleFeature feature) { try { - Polygon p = (Polygon) feature.getDefaultGeometry(); - Point centroid = p.getCentroid(); + Polygon polygon = (Polygon) feature.getDefaultGeometry(); + Point centroid = polygon.getCentroid(); +// CoordinateReferenceSystem crs = feature.getDefaultGeometryProperty().getType() +// .getCoordinateReferenceSystem(); + // TODO convert the centroid +// if (!crs.getName().equals(DefaultGeographicCRS.WGS84.getName())) +// throw new IllegalArgumentException("Only WGS84 CRS is currently supported"); + + // create automatic CRS for optimal computation String code = "AUTO:42001," + centroid.getX() + "," + centroid.getY(); CoordinateReferenceSystem auto = CRS.decode(code); MathTransform transform = CRS.findMathTransform(DefaultGeographicCRS.WGS84, auto); - Polygon projed = (Polygon) JTS.transform(p, transform); + Polygon projed = (Polygon) JTS.transform(polygon, transform); return Quantities.getQuantity(projed.getArea(), SI.SQUARE_METRE); } catch (MismatchedDimensionException | FactoryException | TransformException e) { - throw new IllegalStateException("Cannot claculate area of feature"); + throw new IllegalStateException("Cannot calculate area of feature", e); + } + } + + /** In square meters. The {@link Polygon} must use WGS84 coordinates. */ + public static Quantity calcArea(Polygon polygon) { + assert polygon.getSRID() == 0 || polygon.getSRID() == 4326; + return calcArea(polygon, DefaultGeographicCRS.WGS84, polygon.getCentroid()); + } + + public static Quantity calcArea(Polygon polygon, CoordinateReferenceSystem crs, Point wgs84centroid) { + try { + // create automatic CRS for optimal computation + String code = "AUTO:42001," + wgs84centroid.getX() + "," + wgs84centroid.getY(); + CoordinateReferenceSystem auto = CRS.decode(code); + + MathTransform transform = CRS.findMathTransform(crs, auto); + + Polygon projed = (Polygon) JTS.transform(polygon, transform); + return Quantities.getQuantity(projed.getArea(), SI.SQUARE_METRE); + } catch (MismatchedDimensionException | FactoryException | TransformException e) { + throw new IllegalStateException("Cannot calculate area of feature", e); } } @@ -342,13 +375,81 @@ public class GeoUtils { } public static long getScaleFromResolution(long resolution) { - // see https://gis.stackexchange.com/questions/242424/how-to-get-map-units-to-find-current-scale-in-openlayers + // see + // https://gis.stackexchange.com/questions/242424/how-to-get-map-units-to-find-current-scale-in-openlayers final double INCHES_PER_UNIT = 39.37;// m // final double INCHES_PER_UNIT = 4374754;// dd final long DOTS_PER_INCH = 72; return Math.round(INCHES_PER_UNIT * DOTS_PER_INCH * resolution); } + /* + * GEOMETRY + */ + /** + * Ensure that a {@link Polygon} is valid and simple by removing artefacts + * (typically coming from GPS). + * + * @return a simple and valid polygon, or null if the ignoredArea ratio is above + * the provided threshold. + */ + public static Polygon cleanPolygon(Polygon polygon, double ignoredAreaRatio) { + Polygonizer polygonizer = new Polygonizer(true); + Geometry fixed = GeometryFixer.fix(polygon, false); + polygonizer.add(fixed); + @SuppressWarnings("unchecked") + List polygons = new ArrayList<>(polygonizer.getPolygons()); + + if (polygons.size() == 0) { + throw new IllegalStateException("Polygonizer failed to extract any polygon"); + } + + Polygon best; + if (polygons.size() == 1) { + best = polygons.get(0); + } else { + double totalArea = fixed.getArea(); + best = polygons.get(0); + double bestArea = best.getArea(); + for (int i = 1; i < polygons.size(); i++) { + Polygon p = polygons.get(i); + double a = p.getArea(); + if (a > bestArea) { + best = p; + bestArea = a; + } else { + // double ratio = a / totalArea; + } + } + double ignoredRatio = (totalArea - bestArea) / totalArea; + if (ignoredRatio > ignoredAreaRatio) + return null; + + if (!best.isValid() || !best.isSimple()) { + throw new IllegalStateException("The polygon is not simple and/or valid after cleaning"); + } + } + // while we are here, we make sure that the geometry will be normalised + best.normalize(); + return best; + } + + /** + * The smallest polygon without holes containing all the points in this + * geometry. + */ + public static Polygon concaveHull(Geometry geom, double lengthRatio) { + Objects.requireNonNull(geom); + if (geom.getNumPoints() < 3) + throw new IllegalArgumentException("At least 3 points are reuired to compute the concave hull and geometry " + + geom.getClass() + " contains only " + geom.getNumPoints()); + Geometry hull = ConcaveHull.concaveHullByLengthRatio(geom, lengthRatio, false); + if (hull instanceof Polygon polygon) + return polygon; + else + throw new IllegalStateException("Hull is not a polygon but a " + hull.getClass()); + } + /** Singleton. */ private GeoUtils() { } diff --git a/org.argeo.app.geo/src/org/argeo/app/geo/GpxUtils.java b/org.argeo.app.geo/src/org/argeo/app/geo/GpxUtils.java index 2bcff15..5b07d16 100644 --- a/org.argeo.app.geo/src/org/argeo/app/geo/GpxUtils.java +++ b/org.argeo.app.geo/src/org/argeo/app/geo/GpxUtils.java @@ -22,6 +22,7 @@ import org.locationtech.jts.geom.Coordinate; import org.locationtech.jts.geom.Geometry; import org.locationtech.jts.geom.GeometryFactory; import org.locationtech.jts.geom.LineString; +import org.locationtech.jts.geom.MultiPoint; import org.locationtech.jts.geom.Polygon; import org.opengis.feature.simple.SimpleFeature; import org.opengis.feature.simple.SimpleFeatureType; @@ -52,7 +53,7 @@ public class GpxUtils { */ @SuppressWarnings("unchecked") public static T parseGpxTrackTo(InputStream in, Class clss) throws IOException { - GeometryFactory geometryFactory = JTS.GEOMETRY_FACTORY; + GeometryFactory geometryFactory = JTS.GEOMETRY_FACTORY_WGS84; List coordinates = new ArrayList<>(); try { SAXParserFactory factory = SAXParserFactory.newInstance(); @@ -81,14 +82,17 @@ public class GpxUtils { LineString lineString = geometryFactory .createLineString(coordinates.toArray(new Coordinate[coordinates.size()])); return (T) lineString; + } else if (MultiPoint.class.isAssignableFrom(clss)) { + MultiPoint multiPoint = geometryFactory + .createMultiPointFromCoords(coordinates.toArray(new Coordinate[coordinates.size()])); + // multiPoint.normalize(); + return (T) multiPoint; } else if (Polygon.class.isAssignableFrom(clss)) { // close the line string coordinates.add(coordinates.get(0)); Polygon polygon = geometryFactory.createPolygon(coordinates.toArray(new Coordinate[coordinates.size()])); return (T) polygon; - } - // TODO MultiPoint? MultiLine? etc. - else if (SimpleFeature.class.isAssignableFrom(clss)) { + } else if (SimpleFeature.class.isAssignableFrom(clss)) { SimpleFeatureBuilder featureBuilder = new SimpleFeatureBuilder(LINESTRING_FEATURE_TYPE); LineString lineString = geometryFactory .createLineString(coordinates.toArray(new Coordinate[coordinates.size()])); @@ -99,7 +103,7 @@ public class GpxUtils { throw new IllegalArgumentException("Unsupported format " + clss); } } - + /** @deprecated Use {@link #parseGpxTrackTo(InputStream, Class)} instead. */ @Deprecated public static SimpleFeature parseGpxToPolygon(InputStream in) throws IOException { diff --git a/org.argeo.app.geo/src/org/argeo/app/geo/http/WfsHttpHandler.java b/org.argeo.app.geo/src/org/argeo/app/geo/http/WfsHttpHandler.java index 4ab3654..84b9893 100644 --- a/org.argeo.app.geo/src/org/argeo/app/geo/http/WfsHttpHandler.java +++ b/org.argeo.app.geo/src/org/argeo/app/geo/http/WfsHttpHandler.java @@ -28,7 +28,7 @@ import org.argeo.app.api.EntityType; import org.argeo.app.api.WGS84PosName; import org.argeo.app.api.geo.FeatureAdapter; import org.argeo.app.geo.CqlUtils; -import org.argeo.app.geo.GeoJSon; +import org.argeo.app.geo.GeoJson; import org.argeo.app.geo.GpxUtils; import org.argeo.app.geo.JTS; import org.argeo.cms.acr.json.AcrJsonUtils; @@ -201,12 +201,12 @@ public class WfsHttpHandler implements HttpHandler { String featureId = getFeatureId(c); if (featureId != null) generator.write("id", featureId); - GeoJSon.writeBBox(generator, defaultGeometry); - generator.writeStartObject(GeoJSon.GEOMETRY); - GeoJSon.writeGeometry(generator, defaultGeometry); +// GeoJson.writeBBox(generator, defaultGeometry); + generator.writeStartObject(GeoJson.GEOMETRY); + GeoJson.writeGeometry(generator, defaultGeometry); generator.writeEnd();// geometry object - generator.writeStartObject(GeoJSon.PROPERTIES); + generator.writeStartObject(GeoJson.PROPERTIES); AcrJsonUtils.writeTimeProperties(generator, c); if (featureAdapter != null) featureAdapter.writeProperties(generator, c, typeName); -- 2.30.2