X-Git-Url: http://git.argeo.org/?a=blobdiff_plain;ds=inline;f=org.argeo.app.geo%2Fsrc%2Forg%2Fargeo%2Fapp%2Fgeo%2FGeoUtils.java;h=f13467821de3984e37a61e48b0458f6c97858d54;hb=ee238f9b4c51d5a7e372441d680c75931694569f;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() {
}