X-Git-Url: http://git.argeo.org/?a=blobdiff_plain;f=org.argeo.app.geo%2Fsrc%2Forg%2Fargeo%2Fapp%2Fgeo%2FGeoUtils.java;h=e565866bc123cb89ce9ddc172d84b4ee0bd84f7a;hb=03c171b5d502b461c534816c3c11cae889aef3c1;hp=3cec484c2f14d358f86b2cb66bbcab8bb8ee51c4;hpb=737346afd15e56f9339a7c41ed4e26d65bbcbe69;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 3cec484..e565866 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,76 +8,111 @@ 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.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); + } + } + + /** 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); } } @@ -176,10 +211,14 @@ public class GeoUtils { } } - public static org.opengis.style.Style createStyleFromCss(String css) { + /* + * STYLING + */ + + 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(); @@ -190,31 +229,32 @@ public class GeoUtils { // e.printStackTrace(); // } - return style; + return (Style) style; } public static String createSldFromCss(String name, String title, String css) { - - StyleFactory sf = CommonFactoryFinder.getStyleFactory(); - - StyledLayerDescriptor sld = sf.createStyledLayerDescriptor(); + StyledLayerDescriptor sld = GeoTools.STYLE_FACTORY.createStyledLayerDescriptor(); sld.setName(name); sld.setTitle(title); - UserLayer layer = sf.createUserLayer(); - layer.setName("default"); + NamedLayer layer = GeoTools.STYLE_FACTORY.createNamedLayer(); + layer.setName(name); - org.opengis.style.Style style = createStyleFromCss(css); - layer.userStyles().add((Style) style); + Style style = createStyleFromCss(css); + layer.styles().add(style); sld.layers().add(layer); + return sldToXml(sld); + } + + public static String sldToXml(StyledLayerDescriptor sld) { try { SLDTransformer styleTransform = new SLDTransformer(); String xml = styleTransform.transform(sld); // System.out.println(xml); return xml; } catch (TransformerException e) { - throw new IllegalStateException(e); + throw new IllegalStateException("Cannot write SLD as XML", e); } } @@ -236,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"); @@ -335,6 +375,82 @@ 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 + 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() { }