1 package org
.argeo
.app
.geo
;
3 import java
.io
.IOException
;
4 import java
.io
.Serializable
;
6 import java
.nio
.file
.Path
;
7 import java
.util
.ArrayList
;
8 import java
.util
.HashMap
;
11 import java
.util
.Objects
;
13 import javax
.measure
.Quantity
;
14 import javax
.measure
.quantity
.Area
;
15 import javax
.xml
.transform
.TransformerException
;
17 import org
.geotools
.api
.data
.SimpleFeatureSource
;
18 import org
.geotools
.api
.data
.SimpleFeatureStore
;
19 import org
.geotools
.api
.data
.Transaction
;
20 import org
.geotools
.api
.feature
.simple
.SimpleFeature
;
21 import org
.geotools
.api
.feature
.simple
.SimpleFeatureType
;
22 import org
.geotools
.api
.filter
.Filter
;
23 import org
.geotools
.api
.filter
.FilterFactory
;
24 import org
.geotools
.api
.filter
.expression
.Expression
;
25 import org
.geotools
.api
.geometry
.MismatchedDimensionException
;
26 import org
.geotools
.api
.referencing
.FactoryException
;
27 import org
.geotools
.api
.referencing
.crs
.CoordinateReferenceSystem
;
28 import org
.geotools
.api
.referencing
.operation
.MathTransform
;
29 import org
.geotools
.api
.referencing
.operation
.TransformException
;
30 import org
.geotools
.api
.style
.AnchorPoint
;
31 import org
.geotools
.api
.style
.Displacement
;
32 import org
.geotools
.api
.style
.FeatureTypeConstraint
;
33 import org
.geotools
.api
.style
.FeatureTypeStyle
;
34 import org
.geotools
.api
.style
.Fill
;
35 import org
.geotools
.api
.style
.Graphic
;
36 import org
.geotools
.api
.style
.GraphicalSymbol
;
37 import org
.geotools
.api
.style
.NamedLayer
;
38 import org
.geotools
.api
.style
.PointSymbolizer
;
39 import org
.geotools
.api
.style
.Rule
;
40 import org
.geotools
.api
.style
.Stroke
;
41 import org
.geotools
.api
.style
.Style
;
42 import org
.geotools
.api
.style
.StyleFactory
;
43 import org
.geotools
.api
.style
.StyledLayerDescriptor
;
44 import org
.geotools
.api
.style
.UserLayer
;
45 import org
.geotools
.data
.DefaultTransaction
;
46 import org
.geotools
.data
.collection
.ListFeatureCollection
;
47 import org
.geotools
.data
.shapefile
.ShapefileDataStore
;
48 import org
.geotools
.data
.shapefile
.ShapefileDataStoreFactory
;
49 import org
.geotools
.data
.simple
.SimpleFeatureCollection
;
50 import org
.geotools
.factory
.CommonFactoryFinder
;
51 import org
.geotools
.geometry
.jts
.JTS
;
52 import org
.geotools
.referencing
.CRS
;
53 import org
.geotools
.referencing
.crs
.DefaultGeographicCRS
;
54 import org
.geotools
.styling
.css
.CssParser
;
55 import org
.geotools
.styling
.css
.CssTranslator
;
56 import org
.geotools
.styling
.css
.Stylesheet
;
57 import org
.geotools
.xml
.styling
.SLDTransformer
;
58 import org
.locationtech
.jts
.algorithm
.hull
.ConcaveHull
;
59 import org
.locationtech
.jts
.geom
.Coordinate
;
60 import org
.locationtech
.jts
.geom
.Geometry
;
61 import org
.locationtech
.jts
.geom
.Point
;
62 import org
.locationtech
.jts
.geom
.Polygon
;
63 import org
.locationtech
.jts
.geom
.util
.GeometryFixer
;
64 import org
.locationtech
.jts
.operation
.polygonize
.Polygonizer
;
67 import tech
.units
.indriya
.quantity
.Quantities
;
69 /** Utilities around geographical format, mostly wrapping GeoTools patterns. */
70 public class GeoUtils
{
71 public final static String EPSG_4326
= "EPSG:4326";
73 /** In square meters. */
74 public static Quantity
<Area
> calcArea(SimpleFeature feature
) {
76 Polygon polygon
= (Polygon
) feature
.getDefaultGeometry();
77 Point centroid
= polygon
.getCentroid();
78 // CoordinateReferenceSystem crs = feature.getDefaultGeometryProperty().getType()
79 // .getCoordinateReferenceSystem();
80 // TODO convert the centroid
81 // if (!crs.getName().equals(DefaultGeographicCRS.WGS84.getName()))
82 // throw new IllegalArgumentException("Only WGS84 CRS is currently supported");
84 // create automatic CRS for optimal computation
85 String code
= "AUTO:42001," + centroid
.getX() + "," + centroid
.getY();
86 CoordinateReferenceSystem auto
= CRS
.decode(code
);
88 MathTransform transform
= CRS
.findMathTransform(DefaultGeographicCRS
.WGS84
, auto
);
90 Polygon projed
= (Polygon
) JTS
.transform(polygon
, transform
);
91 return Quantities
.getQuantity(projed
.getArea(), SI
.SQUARE_METRE
);
92 } catch (MismatchedDimensionException
| FactoryException
| TransformException e
) {
93 throw new IllegalStateException("Cannot calculate area of feature", e
);
97 /** In square meters. The {@link Polygon} must use WGS84 coordinates. */
98 public static Quantity
<Area
> calcArea(Polygon polygon
) {
99 assert polygon
.getSRID() == 0 || polygon
.getSRID() == 4326;
100 return calcArea(polygon
, DefaultGeographicCRS
.WGS84
, polygon
.getCentroid());
103 public static Quantity
<Area
> calcArea(Polygon polygon
, CoordinateReferenceSystem crs
, Point wgs84centroid
) {
105 // create automatic CRS for optimal computation
106 String code
= "AUTO:42001," + wgs84centroid
.getX() + "," + wgs84centroid
.getY();
107 CoordinateReferenceSystem auto
= CRS
.decode(code
);
109 MathTransform transform
= CRS
.findMathTransform(crs
, auto
);
111 Polygon projed
= (Polygon
) JTS
.transform(polygon
, transform
);
112 return Quantities
.getQuantity(projed
.getArea(), SI
.SQUARE_METRE
);
113 } catch (MismatchedDimensionException
| FactoryException
| TransformException e
) {
114 throw new IllegalStateException("Cannot calculate area of feature", e
);
118 public static void exportToSvg(Geometry
[] geometries
, Writer out
, int width
, int height
) {
120 double minY
= Double
.POSITIVE_INFINITY
;
121 double maxY
= Double
.NEGATIVE_INFINITY
;
122 double minX
= Double
.POSITIVE_INFINITY
;
123 double maxX
= Double
.NEGATIVE_INFINITY
;
124 List
<String
> shapes
= new ArrayList
<>();
125 for (Geometry geometry
: geometries
) {
126 StringBuffer sb
= new StringBuffer();
127 sb
.append("<polyline style=\"stroke-width:1;stroke:#000000;fill:none;\" points=\"");
129 if (geometry
instanceof Polygon p
) {
130 Point centroid
= p
.getCentroid();
131 String code
= "AUTO:42001," + centroid
.getX() + "," + centroid
.getY();
132 CoordinateReferenceSystem auto
= CRS
.decode(code
);
134 MathTransform transform
= CRS
.findMathTransform(DefaultGeographicCRS
.WGS84
, auto
);
136 Polygon projed
= (Polygon
) JTS
.transform(p
, transform
);
138 // EPSG4326 are in lat/lon order, so we invert coordinates
139 for (Coordinate coord
: projed
.getCoordinates()) {
150 sb
.append(x
+ "," + y
+ " ");
153 sb
.append("</polyline>\n");
154 shapes
.add(sb
.toString());
156 throw new IllegalArgumentException("Unsuppported geometry type " + geometry
.getClass().getName());
159 double viewportHeight
= maxY
- minY
;
160 double viewportWidth
= maxX
- minX
;
161 out
.write("<svg xmlns=\"http://www.w3.org/2000/svg\"\n");
162 out
.write(" width=\"" + width
+ "\"\n");
163 out
.write(" height=\"" + height
+ "\"\n");
164 out
.write(" viewBox=\"" + minX
+ " " + minY
+ " " + viewportWidth
+ " " + viewportHeight
+ "\"\n");
165 out
.write(" preserveAspectRatio=\"xMidYMid meet\"\n");
167 for (String shape
: shapes
) {
172 } catch (IOException
| FactoryException
| MismatchedDimensionException
| TransformException e
) {
173 throw new RuntimeException("Cannot export to SVG", e
);
177 /** Write a list of simple features to a shapefile. */
178 public static void saveFeaturesAsShapefile(SimpleFeatureType featureType
, List
<SimpleFeature
> features
,
181 ShapefileDataStoreFactory dataStoreFactory
= new ShapefileDataStoreFactory();
183 Map
<String
, Serializable
> params
= new HashMap
<>();
184 params
.put("url", shpFile
.toUri().toURL());
186 params
.put("create spatial index", Boolean
.TRUE
);
188 ShapefileDataStore newDataStore
= (ShapefileDataStore
) dataStoreFactory
.createNewDataStore(params
);
189 newDataStore
.createSchema(featureType
);
191 String typeName
= newDataStore
.getTypeNames()[0];
192 SimpleFeatureSource featureSource
= newDataStore
.getFeatureSource(typeName
);
193 if (featureSource
instanceof SimpleFeatureStore
) {
194 SimpleFeatureStore featureStore
= (SimpleFeatureStore
) featureSource
;
195 SimpleFeatureCollection collection
= new ListFeatureCollection(featureType
, features
);
197 try (Transaction transaction
= new DefaultTransaction("create")) {
199 featureStore
.setTransaction(transaction
);
200 featureStore
.addFeatures(collection
);
201 transaction
.commit();
202 } catch (Exception problem
) {
203 transaction
.rollback();
204 throw new RuntimeException("Cannot write shapefile " + shpFile
, problem
);
208 throw new IllegalArgumentException(typeName
+ " does not support read/write access");
210 } catch (IOException e
) {
211 throw new RuntimeException("Cannot write shapefile " + shpFile
, e
);
219 public static Style
createStyleFromCss(String css
) {
220 Stylesheet ss
= CssParser
.parse(css
);
221 CssTranslator translator
= new CssTranslator();
222 Style style
= translator
.translate(ss
);
225 // SLDTransformer styleTransform = new SLDTransformer();
226 // String xml = styleTransform.transform(style);
227 // System.out.println(xml);
228 // } catch (TransformerException e) {
229 // // TODO Auto-generated catch block
230 // e.printStackTrace();
233 return (Style
) style
;
236 public static String
createSldFromCss(String name
, String title
, String css
) {
237 StyledLayerDescriptor sld
= GeoTools
.STYLE_FACTORY
.createStyledLayerDescriptor();
241 NamedLayer layer
= GeoTools
.STYLE_FACTORY
.createNamedLayer();
244 Style style
= createStyleFromCss(css
);
245 layer
.styles().add(style
);
247 sld
.layers().add(layer
);
248 return sldToXml(sld
);
251 public static String
sldToXml(StyledLayerDescriptor sld
) {
253 SLDTransformer styleTransform
= new SLDTransformer();
254 String xml
= styleTransform
.transform(sld
);
255 // System.out.println(xml);
257 } catch (TransformerException e
) {
258 throw new IllegalStateException("Cannot write SLD as XML", e
);
262 public static void main(String
... args
) {
265 mark: symbol(circle);
274 createSldFromCss("test", "Test", css
);
277 public static String
createTestSLD() {
279 StyleFactory sf
= CommonFactoryFinder
.getStyleFactory();
280 FilterFactory ff
= CommonFactoryFinder
.getFilterFactory();
282 StyledLayerDescriptor sld
= sf
.createStyledLayerDescriptor();
284 sld
.setTitle("Example");
285 sld
.setAbstract("Example Style Layer Descriptor");
287 UserLayer layer
= sf
.createUserLayer();
288 layer
.setName("layer");
291 // define constraint limited what features the sld applies to
292 FeatureTypeConstraint constraint
= sf
.createFeatureTypeConstraint("Feature", Filter
.INCLUDE
);
294 layer
.layerFeatureConstraints().add(constraint
);
297 // create a "user defined" style
298 Style style
= sf
.createStyle();
299 style
.setName("style");
300 style
.getDescription().setTitle("User Style");
301 style
.getDescription().setAbstract("Definition of Style");
304 // define feature type styles used to actually define how features are rendered
305 FeatureTypeStyle featureTypeStyle
= sf
.createFeatureTypeStyle();
308 // first rule to draw cities
309 Rule rule1
= sf
.createRule();
310 rule1
.setName("rule1");
311 rule1
.getDescription().setTitle("City");
312 rule1
.getDescription().setAbstract("Rule for drawing cities");
313 // rule1.setFilter(ff.less(ff.property("POPULATION"), ff.literal(50000)));
316 // create the graphical mark used to represent a city
317 Stroke stroke
= sf
.stroke(ff
.literal("#000000"), null, null, null, null, null, null);
318 Fill fill
= sf
.fill(null, ff
.literal(java
.awt
.Color
.BLUE
), ff
.literal(1.0));
320 // // OnLineResource implemented by gt-metadata - so no factory!
321 // OnLineResourceImpl svg = new OnLineResourceImpl(new URI("file:city.svg"));
322 // svg.freeze(); // freeze to prevent modification at runtime
324 // OnLineResourceImpl png = new OnLineResourceImpl(new URI("file:city.png"));
325 // png.freeze(); // freeze to prevent modification at runtime
328 // List of symbols is considered in order with the rendering engine choosing
329 // the first one it can handle. Allowing for svg, png, mark order
330 List
<GraphicalSymbol
> symbols
= new ArrayList
<>();
331 // symbols.add(sf.externalGraphic(svg, "svg", null)); // svg preferred
332 // symbols.add(sf.externalGraphic(png, "png", null)); // png preferred
333 symbols
.add(sf
.mark(ff
.literal("circle"), fill
, stroke
)); // simple circle backup plan
335 Expression opacity
= null; // use default
336 Expression size
= ff
.literal(10);
337 Expression rotation
= null; // use default
338 AnchorPoint anchor
= null; // use default
339 Displacement displacement
= null; // use default
341 // define a point symbolizer of a small circle
342 Graphic city
= sf
.graphic(symbols
, opacity
, size
, rotation
, anchor
, displacement
);
343 PointSymbolizer pointSymbolizer
= sf
.pointSymbolizer("point", ff
.property("the_geom"), null, null, city
);
345 rule1
.symbolizers().add(pointSymbolizer
);
347 featureTypeStyle
.rules().add(rule1
);
352 // List<GraphicalSymbol> dotSymbols = new ArrayList<>();
353 // dotSymbols.add(sf.mark(ff.literal("circle"), null, null));
354 // Graphic dotGraphic = sf.graphic(dotSymbols, null, ff.literal(3), null, null, null);
355 // PointSymbolizer dotSymbolizer = sf.pointSymbolizer("dot", ff.property("the_geom"), null, null, dotGraphic);
356 // List<org.opengis.style.Symbolizer> symbolizers = new ArrayList<>();
357 // symbolizers.add(dotSymbolizer);
358 // Filter other = null; // null will mark this rule as "other" accepting all remaining features
359 // Rule rule2 = sf.rule("default", null, null, Double.MIN_VALUE, Double.MAX_VALUE, symbolizers, other);
360 // featureTypeStyle.rules().add(rule2);
362 style
.featureTypeStyles().add(featureTypeStyle
);
364 layer
.userStyles().add(style
);
366 sld
.layers().add(layer
);
369 SLDTransformer styleTransform
= new SLDTransformer();
370 String xml
= styleTransform
.transform(sld
);
371 System
.out
.println(xml
);
373 } catch (TransformerException e
) {
374 throw new IllegalStateException(e
);
379 public static long getScaleFromResolution(long resolution
) {
381 // https://gis.stackexchange.com/questions/242424/how-to-get-map-units-to-find-current-scale-in-openlayers
382 final double INCHES_PER_UNIT
= 39.37;// m
383 // final double INCHES_PER_UNIT = 4374754;// dd
384 final long DOTS_PER_INCH
= 72;
385 return Math
.round(INCHES_PER_UNIT
* DOTS_PER_INCH
* resolution
);
392 * Ensure that a {@link Polygon} is valid and simple by removing artefacts
393 * (typically coming from GPS).
395 * @return a simple and valid polygon, or null if the ignoredArea ratio is above
396 * the provided threshold.
398 public static Polygon
cleanPolygon(Polygon polygon
, double ignoredAreaRatio
) {
399 Polygonizer polygonizer
= new Polygonizer(true);
400 Geometry fixed
= GeometryFixer
.fix(polygon
, false);
401 polygonizer
.add(fixed
);
402 @SuppressWarnings("unchecked")
403 List
<Polygon
> polygons
= new ArrayList
<>(polygonizer
.getPolygons());
405 if (polygons
.size() == 0) {
406 throw new IllegalStateException("Polygonizer failed to extract any polygon");
410 if (polygons
.size() == 1) {
411 best
= polygons
.get(0);
413 double totalArea
= fixed
.getArea();
414 best
= polygons
.get(0);
415 double bestArea
= best
.getArea();
416 for (int i
= 1; i
< polygons
.size(); i
++) {
417 Polygon p
= polygons
.get(i
);
418 double a
= p
.getArea();
423 // double ratio = a / totalArea;
426 double ignoredRatio
= (totalArea
- bestArea
) / totalArea
;
427 if (ignoredRatio
> ignoredAreaRatio
)
430 if (!best
.isValid() || !best
.isSimple()) {
431 throw new IllegalStateException("The polygon is not simple and/or valid after cleaning");
434 // while we are here, we make sure that the geometry will be normalised
440 * The smallest polygon without holes containing all the points in this
443 public static Polygon
concaveHull(Geometry geom
, double lengthRatio
) {
444 Objects
.requireNonNull(geom
);
445 if (geom
.getNumPoints() < 3)
446 throw new IllegalArgumentException("At least 3 points are reuired to compute the concave hull and geometry "
447 + geom
.getClass() + " contains only " + geom
.getNumPoints());
448 Geometry hull
= ConcaveHull
.concaveHullByLengthRatio(geom
, lengthRatio
, false);
449 if (hull
instanceof Polygon polygon
)
452 throw new IllegalStateException("Hull is not a polygon but a " + hull
.getClass());