Introduce geometry to SVG
[gpl/argeo-suite.git] / org.argeo.app.geo / src / org / argeo / app / geo / GeoUtils.java
index 8f971291aeb49c877e145bee691ff92ea87dac51..47096ccc0d1f5b7b2fddb2ce27cedf9d705235b6 100644 (file)
@@ -8,118 +8,152 @@ 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;
 import javax.xml.transform.TransformerException;
 
+import org.geotools.api.data.SimpleFeatureSource;
+import org.geotools.api.data.SimpleFeatureStore;
+import org.geotools.api.data.Transaction;
+import org.geotools.api.feature.simple.SimpleFeature;
+import org.geotools.api.feature.simple.SimpleFeatureType;
+import org.geotools.api.filter.Filter;
+import org.geotools.api.filter.FilterFactory;
+import org.geotools.api.filter.expression.Expression;
+import org.geotools.api.geometry.MismatchedDimensionException;
+import org.geotools.api.referencing.FactoryException;
+import org.geotools.api.referencing.crs.CoordinateReferenceSystem;
+import org.geotools.api.referencing.operation.MathTransform;
+import org.geotools.api.referencing.operation.TransformException;
+import org.geotools.api.style.AnchorPoint;
+import org.geotools.api.style.Displacement;
+import org.geotools.api.style.FeatureTypeConstraint;
+import org.geotools.api.style.FeatureTypeStyle;
+import org.geotools.api.style.Fill;
+import org.geotools.api.style.Graphic;
+import org.geotools.api.style.GraphicalSymbol;
+import org.geotools.api.style.NamedLayer;
+import org.geotools.api.style.PointSymbolizer;
+import org.geotools.api.style.Rule;
+import org.geotools.api.style.Stroke;
+import org.geotools.api.style.Style;
+import org.geotools.api.style.StyleFactory;
+import org.geotools.api.style.StyledLayerDescriptor;
+import org.geotools.api.style.UserLayer;
 import org.geotools.data.DefaultTransaction;
-import org.geotools.data.Transaction;
 import org.geotools.data.collection.ListFeatureCollection;
 import org.geotools.data.shapefile.ShapefileDataStore;
 import org.geotools.data.shapefile.ShapefileDataStoreFactory;
 import org.geotools.data.simple.SimpleFeatureCollection;
-import org.geotools.data.simple.SimpleFeatureIterator;
-import org.geotools.data.simple.SimpleFeatureSource;
-import org.geotools.data.simple.SimpleFeatureStore;
 import org.geotools.factory.CommonFactoryFinder;
 import org.geotools.geometry.jts.JTS;
 import org.geotools.referencing.CRS;
 import org.geotools.referencing.crs.DefaultGeographicCRS;
-import org.geotools.styling.AnchorPoint;
-import org.geotools.styling.Displacement;
-import org.geotools.styling.FeatureTypeConstraint;
-import org.geotools.styling.FeatureTypeStyle;
-import org.geotools.styling.Fill;
-import org.geotools.styling.Graphic;
-import org.geotools.styling.NamedLayer;
-import org.geotools.styling.PointSymbolizer;
-import org.geotools.styling.Rule;
-import org.geotools.styling.Stroke;
-import org.geotools.styling.Style;
-import org.geotools.styling.StyleFactory;
-import org.geotools.styling.StyledLayerDescriptor;
-import org.geotools.styling.UserLayer;
 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.opengis.feature.simple.SimpleFeature;
-import org.opengis.feature.simple.SimpleFeatureType;
-import org.opengis.filter.Filter;
-import org.opengis.filter.FilterFactory2;
-import org.opengis.filter.expression.Expression;
-import org.opengis.geometry.MismatchedDimensionException;
-import org.opengis.referencing.FactoryException;
-import org.opengis.referencing.crs.CoordinateReferenceSystem;
-import org.opengis.referencing.operation.MathTransform;
-import org.opengis.referencing.operation.TransformException;
-import org.opengis.style.GraphicalSymbol;
+import org.locationtech.jts.geom.util.GeometryFixer;
+import org.locationtech.jts.operation.polygonize.Polygonizer;
 
 import si.uom.SI;
 import tech.units.indriya.quantity.Quantities;
 
 /** Utilities around geographical format, mostly wrapping GeoTools patterns. */
 public class GeoUtils {
+       public final static String EPSG_4326 = "EPSG:4326";
 
        /** In square meters. */
        public static Quantity<Area> 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);
                }
        }
 
-       public static void exportToSvg(SimpleFeatureCollection features, Writer out, int width, int height) {
+       /** In square meters. The {@link Polygon} must use WGS84 coordinates. */
+       public static Quantity<Area> calcArea(Polygon polygon) {
+               assert polygon.getSRID() == 0 || polygon.getSRID() == 4326;
+               return calcArea(polygon, DefaultGeographicCRS.WGS84, polygon.getCentroid());
+       }
+
+       public static Quantity<Area> 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);
+               }
+       }
+
+       public static void exportToSvg(Geometry[] geometries, Writer out, int width, int height) {
                try {
                        double minY = Double.POSITIVE_INFINITY;
                        double maxY = Double.NEGATIVE_INFINITY;
                        double minX = Double.POSITIVE_INFINITY;
                        double maxX = Double.NEGATIVE_INFINITY;
                        List<String> shapes = new ArrayList<>();
-                       for (SimpleFeatureIterator it = features.features(); it.hasNext();) {
-                               SimpleFeature feature = it.next();
+                       for (Geometry geometry : geometries) {
                                StringBuffer sb = new StringBuffer();
                                sb.append("<polyline style=\"stroke-width:1;stroke:#000000;fill:none;\" points=\"");
 
-                               Polygon p = (Polygon) feature.getDefaultGeometry();
-                               Point centroid = p.getCentroid();
-                               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);
-
-                               for (Coordinate coord : projed.getCoordinates()) {
-                                       double x = coord.x;
-                                       if (x < minX)
-                                               minX = x;
-                                       if (x > maxX)
-                                               maxX = x;
-                                       double y = -coord.y;
-                                       if (y < minY)
-                                               minY = y;
-                                       if (y > maxY)
-                                               maxY = y;
-                                       sb.append(x + "," + y + " ");
+                               if (geometry instanceof Polygon p) {
+                                       Point centroid = p.getCentroid();
+                                       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);
+
+                                       for (Coordinate coord : projed.getCoordinates()) {
+                                               double x = coord.x;
+                                               if (x < minX)
+                                                       minX = x;
+                                               if (x > maxX)
+                                                       maxX = x;
+                                               double y = -coord.y;
+                                               if (y < minY)
+                                                       minY = y;
+                                               if (y > maxY)
+                                                       maxY = y;
+                                               sb.append(x + "," + y + " ");
+                                       }
+                                       sb.append("\">");
+                                       sb.append("</polyline>\n");
+                                       shapes.add(sb.toString());
+                               } else {
+                                       throw new IllegalArgumentException("Unsuppported geometry type " + geometry.getClass().getName());
                                }
-                               sb.append("\">");
-                               sb.append("</polyline>\n");
-                               shapes.add(sb.toString());
-
                        }
                        double viewportHeight = maxY - minY;
                        double viewportWidth = maxX - minX;
@@ -184,7 +218,7 @@ public class GeoUtils {
        public static Style createStyleFromCss(String css) {
                Stylesheet ss = CssParser.parse(css);
                CssTranslator translator = new CssTranslator();
-               org.opengis.style.Style style = translator.translate(ss);
+               Style style = translator.translate(ss);
 
 //             try {
 //                     SLDTransformer styleTransform = new SLDTransformer();
@@ -242,7 +276,7 @@ public class GeoUtils {
        public static String createTestSLD() {
 
                StyleFactory sf = CommonFactoryFinder.getStyleFactory();
-               FilterFactory2 ff = CommonFactoryFinder.getFilterFactory2();
+               FilterFactory ff = CommonFactoryFinder.getFilterFactory();
 
                StyledLayerDescriptor sld = sf.createStyledLayerDescriptor();
                sld.setName("sld");
@@ -342,13 +376,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<Polygon> 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() {
        }