X-Git-Url: https://git.argeo.org/?a=blobdiff_plain;f=org.argeo.app.geo%2Fsrc%2Forg%2Fargeo%2Fapp%2Fgeo%2FGeoUtils.java;h=f13467821de3984e37a61e48b0458f6c97858d54;hb=2a381ac1f4dddd3d92d105b3684015f526d37f8a;hp=8f971291aeb49c877e145bee691ff92ea87dac51;hpb=8a490e540ac623b3545b1bd3da65ecbf2e4b6436;p=gpl%2Fargeo-suite.git 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..f134678 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,118 +8,153 @@ 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 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 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); + } + } + + 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 shapes = new ArrayList<>(); - for (SimpleFeatureIterator it = features.features(); it.hasNext();) { - SimpleFeature feature = it.next(); + for (Geometry geometry : geometries) { StringBuffer sb = new StringBuffer(); sb.append(" 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); + + // EPSG4326 are in lat/lon order, so we invert coordinates + for (Coordinate coord : projed.getCoordinates()) { + double x = coord.y; + if (x < minX) + minX = x; + if (x > maxX) + maxX = x; + double y = -coord.x; + if (y < minY) + minY = y; + if (y > maxY) + maxY = y; + sb.append(x + "," + y + " "); + } + sb.append("\">"); + sb.append("\n"); + shapes.add(sb.toString()); + } else { + throw new IllegalArgumentException("Unsuppported geometry type " + geometry.getClass().getName()); } - sb.append("\">"); - sb.append("\n"); - shapes.add(sb.toString()); - } double viewportHeight = maxY - minY; double viewportWidth = maxX - minX; @@ -184,7 +219,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 +277,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 +377,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() { }