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
.data
.DefaultTransaction
;
18 import org
.geotools
.data
.Transaction
;
19 import org
.geotools
.data
.collection
.ListFeatureCollection
;
20 import org
.geotools
.data
.shapefile
.ShapefileDataStore
;
21 import org
.geotools
.data
.shapefile
.ShapefileDataStoreFactory
;
22 import org
.geotools
.data
.simple
.SimpleFeatureCollection
;
23 import org
.geotools
.data
.simple
.SimpleFeatureIterator
;
24 import org
.geotools
.data
.simple
.SimpleFeatureSource
;
25 import org
.geotools
.data
.simple
.SimpleFeatureStore
;
26 import org
.geotools
.factory
.CommonFactoryFinder
;
27 import org
.geotools
.geometry
.jts
.JTS
;
28 import org
.geotools
.referencing
.CRS
;
29 import org
.geotools
.referencing
.crs
.DefaultGeographicCRS
;
30 import org
.geotools
.styling
.AnchorPoint
;
31 import org
.geotools
.styling
.Displacement
;
32 import org
.geotools
.styling
.FeatureTypeConstraint
;
33 import org
.geotools
.styling
.FeatureTypeStyle
;
34 import org
.geotools
.styling
.Fill
;
35 import org
.geotools
.styling
.Graphic
;
36 import org
.geotools
.styling
.NamedLayer
;
37 import org
.geotools
.styling
.PointSymbolizer
;
38 import org
.geotools
.styling
.Rule
;
39 import org
.geotools
.styling
.Stroke
;
40 import org
.geotools
.styling
.Style
;
41 import org
.geotools
.styling
.StyleFactory
;
42 import org
.geotools
.styling
.StyledLayerDescriptor
;
43 import org
.geotools
.styling
.UserLayer
;
44 import org
.geotools
.styling
.css
.CssParser
;
45 import org
.geotools
.styling
.css
.CssTranslator
;
46 import org
.geotools
.styling
.css
.Stylesheet
;
47 import org
.geotools
.xml
.styling
.SLDTransformer
;
48 import org
.locationtech
.jts
.algorithm
.hull
.ConcaveHull
;
49 import org
.locationtech
.jts
.geom
.Coordinate
;
50 import org
.locationtech
.jts
.geom
.Geometry
;
51 import org
.locationtech
.jts
.geom
.Point
;
52 import org
.locationtech
.jts
.geom
.Polygon
;
53 import org
.locationtech
.jts
.geom
.util
.GeometryFixer
;
54 import org
.locationtech
.jts
.operation
.polygonize
.Polygonizer
;
55 import org
.opengis
.feature
.simple
.SimpleFeature
;
56 import org
.opengis
.feature
.simple
.SimpleFeatureType
;
57 import org
.opengis
.filter
.Filter
;
58 import org
.opengis
.filter
.FilterFactory2
;
59 import org
.opengis
.filter
.expression
.Expression
;
60 import org
.opengis
.geometry
.MismatchedDimensionException
;
61 import org
.opengis
.referencing
.FactoryException
;
62 import org
.opengis
.referencing
.crs
.CoordinateReferenceSystem
;
63 import org
.opengis
.referencing
.operation
.MathTransform
;
64 import org
.opengis
.referencing
.operation
.TransformException
;
65 import org
.opengis
.style
.GraphicalSymbol
;
68 import tech
.units
.indriya
.quantity
.Quantities
;
70 /** Utilities around geographical format, mostly wrapping GeoTools patterns. */
71 public class GeoUtils
{
72 public final static String EPSG_4326
= "EPSG:4326";
74 /** In square meters. */
75 public static Quantity
<Area
> calcArea(SimpleFeature feature
) {
77 Polygon polygon
= (Polygon
) feature
.getDefaultGeometry();
78 Point centroid
= polygon
.getCentroid();
79 // CoordinateReferenceSystem crs = feature.getDefaultGeometryProperty().getType()
80 // .getCoordinateReferenceSystem();
81 // TODO convert the centroid
82 // if (!crs.getName().equals(DefaultGeographicCRS.WGS84.getName()))
83 // throw new IllegalArgumentException("Only WGS84 CRS is currently supported");
85 // create automatic CRS for optimal computation
86 String code
= "AUTO:42001," + centroid
.getX() + "," + centroid
.getY();
87 CoordinateReferenceSystem auto
= CRS
.decode(code
);
89 MathTransform transform
= CRS
.findMathTransform(DefaultGeographicCRS
.WGS84
, auto
);
91 Polygon projed
= (Polygon
) JTS
.transform(polygon
, transform
);
92 return Quantities
.getQuantity(projed
.getArea(), SI
.SQUARE_METRE
);
93 } catch (MismatchedDimensionException
| FactoryException
| TransformException e
) {
94 throw new IllegalStateException("Cannot calculate area of feature", e
);
98 /** In square meters. The {@link Polygon} must use WGS84 coordinates. */
99 public static Quantity
<Area
> calcArea(Polygon polygon
) {
100 assert polygon
.getSRID() == 0 || polygon
.getSRID() == 4326;
101 return calcArea(polygon
, DefaultGeographicCRS
.WGS84
, polygon
.getCentroid());
104 public static Quantity
<Area
> calcArea(Polygon polygon
, CoordinateReferenceSystem crs
, Point wgs84centroid
) {
106 // create automatic CRS for optimal computation
107 String code
= "AUTO:42001," + wgs84centroid
.getX() + "," + wgs84centroid
.getY();
108 CoordinateReferenceSystem auto
= CRS
.decode(code
);
110 MathTransform transform
= CRS
.findMathTransform(crs
, auto
);
112 Polygon projed
= (Polygon
) JTS
.transform(polygon
, transform
);
113 return Quantities
.getQuantity(projed
.getArea(), SI
.SQUARE_METRE
);
114 } catch (MismatchedDimensionException
| FactoryException
| TransformException e
) {
115 throw new IllegalStateException("Cannot calculate area of feature", e
);
119 public static void exportToSvg(SimpleFeatureCollection features
, Writer out
, int width
, int height
) {
121 double minY
= Double
.POSITIVE_INFINITY
;
122 double maxY
= Double
.NEGATIVE_INFINITY
;
123 double minX
= Double
.POSITIVE_INFINITY
;
124 double maxX
= Double
.NEGATIVE_INFINITY
;
125 List
<String
> shapes
= new ArrayList
<>();
126 for (SimpleFeatureIterator it
= features
.features(); it
.hasNext();) {
127 SimpleFeature feature
= it
.next();
128 StringBuffer sb
= new StringBuffer();
129 sb
.append("<polyline style=\"stroke-width:1;stroke:#000000;fill:none;\" points=\"");
131 Polygon p
= (Polygon
) feature
.getDefaultGeometry();
132 Point centroid
= p
.getCentroid();
133 String code
= "AUTO:42001," + centroid
.getX() + "," + centroid
.getY();
134 CoordinateReferenceSystem auto
= CRS
.decode(code
);
136 MathTransform transform
= CRS
.findMathTransform(DefaultGeographicCRS
.WGS84
, auto
);
138 Polygon projed
= (Polygon
) JTS
.transform(p
, transform
);
140 for (Coordinate coord
: projed
.getCoordinates()) {
151 sb
.append(x
+ "," + y
+ " ");
154 sb
.append("</polyline>\n");
155 shapes
.add(sb
.toString());
158 double viewportHeight
= maxY
- minY
;
159 double viewportWidth
= maxX
- minX
;
160 out
.write("<svg xmlns=\"http://www.w3.org/2000/svg\"\n");
161 out
.write(" width=\"" + width
+ "\"\n");
162 out
.write(" height=\"" + height
+ "\"\n");
163 out
.write(" viewBox=\"" + minX
+ " " + minY
+ " " + viewportWidth
+ " " + viewportHeight
+ "\"\n");
164 out
.write(" preserveAspectRatio=\"xMidYMid meet\"\n");
166 for (String shape
: shapes
) {
171 } catch (IOException
| FactoryException
| MismatchedDimensionException
| TransformException e
) {
172 throw new RuntimeException("Cannot export to SVG", e
);
176 /** Write a list of simple features to a shapefile. */
177 public static void saveFeaturesAsShapefile(SimpleFeatureType featureType
, List
<SimpleFeature
> features
,
180 ShapefileDataStoreFactory dataStoreFactory
= new ShapefileDataStoreFactory();
182 Map
<String
, Serializable
> params
= new HashMap
<>();
183 params
.put("url", shpFile
.toUri().toURL());
185 params
.put("create spatial index", Boolean
.TRUE
);
187 ShapefileDataStore newDataStore
= (ShapefileDataStore
) dataStoreFactory
.createNewDataStore(params
);
188 newDataStore
.createSchema(featureType
);
190 String typeName
= newDataStore
.getTypeNames()[0];
191 SimpleFeatureSource featureSource
= newDataStore
.getFeatureSource(typeName
);
192 if (featureSource
instanceof SimpleFeatureStore
) {
193 SimpleFeatureStore featureStore
= (SimpleFeatureStore
) featureSource
;
194 SimpleFeatureCollection collection
= new ListFeatureCollection(featureType
, features
);
196 try (Transaction transaction
= new DefaultTransaction("create")) {
198 featureStore
.setTransaction(transaction
);
199 featureStore
.addFeatures(collection
);
200 transaction
.commit();
201 } catch (Exception problem
) {
202 transaction
.rollback();
203 throw new RuntimeException("Cannot write shapefile " + shpFile
, problem
);
207 throw new IllegalArgumentException(typeName
+ " does not support read/write access");
209 } catch (IOException e
) {
210 throw new RuntimeException("Cannot write shapefile " + shpFile
, e
);
218 public static Style
createStyleFromCss(String css
) {
219 Stylesheet ss
= CssParser
.parse(css
);
220 CssTranslator translator
= new CssTranslator();
221 org
.opengis
.style
.Style style
= translator
.translate(ss
);
224 // SLDTransformer styleTransform = new SLDTransformer();
225 // String xml = styleTransform.transform(style);
226 // System.out.println(xml);
227 // } catch (TransformerException e) {
228 // // TODO Auto-generated catch block
229 // e.printStackTrace();
232 return (Style
) style
;
235 public static String
createSldFromCss(String name
, String title
, String css
) {
236 StyledLayerDescriptor sld
= GeoTools
.STYLE_FACTORY
.createStyledLayerDescriptor();
240 NamedLayer layer
= GeoTools
.STYLE_FACTORY
.createNamedLayer();
243 Style style
= createStyleFromCss(css
);
244 layer
.styles().add(style
);
246 sld
.layers().add(layer
);
247 return sldToXml(sld
);
250 public static String
sldToXml(StyledLayerDescriptor sld
) {
252 SLDTransformer styleTransform
= new SLDTransformer();
253 String xml
= styleTransform
.transform(sld
);
254 // System.out.println(xml);
256 } catch (TransformerException e
) {
257 throw new IllegalStateException("Cannot write SLD as XML", e
);
261 public static void main(String
... args
) {
264 mark: symbol(circle);
273 createSldFromCss("test", "Test", css
);
276 public static String
createTestSLD() {
278 StyleFactory sf
= CommonFactoryFinder
.getStyleFactory();
279 FilterFactory2 ff
= CommonFactoryFinder
.getFilterFactory2();
281 StyledLayerDescriptor sld
= sf
.createStyledLayerDescriptor();
283 sld
.setTitle("Example");
284 sld
.setAbstract("Example Style Layer Descriptor");
286 UserLayer layer
= sf
.createUserLayer();
287 layer
.setName("layer");
290 // define constraint limited what features the sld applies to
291 FeatureTypeConstraint constraint
= sf
.createFeatureTypeConstraint("Feature", Filter
.INCLUDE
);
293 layer
.layerFeatureConstraints().add(constraint
);
296 // create a "user defined" style
297 Style style
= sf
.createStyle();
298 style
.setName("style");
299 style
.getDescription().setTitle("User Style");
300 style
.getDescription().setAbstract("Definition of Style");
303 // define feature type styles used to actually define how features are rendered
304 FeatureTypeStyle featureTypeStyle
= sf
.createFeatureTypeStyle();
307 // first rule to draw cities
308 Rule rule1
= sf
.createRule();
309 rule1
.setName("rule1");
310 rule1
.getDescription().setTitle("City");
311 rule1
.getDescription().setAbstract("Rule for drawing cities");
312 // rule1.setFilter(ff.less(ff.property("POPULATION"), ff.literal(50000)));
315 // create the graphical mark used to represent a city
316 Stroke stroke
= sf
.stroke(ff
.literal("#000000"), null, null, null, null, null, null);
317 Fill fill
= sf
.fill(null, ff
.literal(java
.awt
.Color
.BLUE
), ff
.literal(1.0));
319 // // OnLineResource implemented by gt-metadata - so no factory!
320 // OnLineResourceImpl svg = new OnLineResourceImpl(new URI("file:city.svg"));
321 // svg.freeze(); // freeze to prevent modification at runtime
323 // OnLineResourceImpl png = new OnLineResourceImpl(new URI("file:city.png"));
324 // png.freeze(); // freeze to prevent modification at runtime
327 // List of symbols is considered in order with the rendering engine choosing
328 // the first one it can handle. Allowing for svg, png, mark order
329 List
<GraphicalSymbol
> symbols
= new ArrayList
<>();
330 // symbols.add(sf.externalGraphic(svg, "svg", null)); // svg preferred
331 // symbols.add(sf.externalGraphic(png, "png", null)); // png preferred
332 symbols
.add(sf
.mark(ff
.literal("circle"), fill
, stroke
)); // simple circle backup plan
334 Expression opacity
= null; // use default
335 Expression size
= ff
.literal(10);
336 Expression rotation
= null; // use default
337 AnchorPoint anchor
= null; // use default
338 Displacement displacement
= null; // use default
340 // define a point symbolizer of a small circle
341 Graphic city
= sf
.graphic(symbols
, opacity
, size
, rotation
, anchor
, displacement
);
342 PointSymbolizer pointSymbolizer
= sf
.pointSymbolizer("point", ff
.property("the_geom"), null, null, city
);
344 rule1
.symbolizers().add(pointSymbolizer
);
346 featureTypeStyle
.rules().add(rule1
);
351 // List<GraphicalSymbol> dotSymbols = new ArrayList<>();
352 // dotSymbols.add(sf.mark(ff.literal("circle"), null, null));
353 // Graphic dotGraphic = sf.graphic(dotSymbols, null, ff.literal(3), null, null, null);
354 // PointSymbolizer dotSymbolizer = sf.pointSymbolizer("dot", ff.property("the_geom"), null, null, dotGraphic);
355 // List<org.opengis.style.Symbolizer> symbolizers = new ArrayList<>();
356 // symbolizers.add(dotSymbolizer);
357 // Filter other = null; // null will mark this rule as "other" accepting all remaining features
358 // Rule rule2 = sf.rule("default", null, null, Double.MIN_VALUE, Double.MAX_VALUE, symbolizers, other);
359 // featureTypeStyle.rules().add(rule2);
361 style
.featureTypeStyles().add(featureTypeStyle
);
363 layer
.userStyles().add(style
);
365 sld
.layers().add(layer
);
368 SLDTransformer styleTransform
= new SLDTransformer();
369 String xml
= styleTransform
.transform(sld
);
370 System
.out
.println(xml
);
372 } catch (TransformerException e
) {
373 throw new IllegalStateException(e
);
378 public static long getScaleFromResolution(long resolution
) {
380 // https://gis.stackexchange.com/questions/242424/how-to-get-map-units-to-find-current-scale-in-openlayers
381 final double INCHES_PER_UNIT
= 39.37;// m
382 // final double INCHES_PER_UNIT = 4374754;// dd
383 final long DOTS_PER_INCH
= 72;
384 return Math
.round(INCHES_PER_UNIT
* DOTS_PER_INCH
* resolution
);
391 * Ensure that a {@link Polygon} is valid and simple by removing artefacts
392 * (typically coming from GPS).
394 * @return a simple and valid polygon, or null if the ignoredArea ratio is above
395 * the provided threshold.
397 public static Polygon
cleanPolygon(Polygon polygon
, double ignoredAreaRatio
) {
398 Polygonizer polygonizer
= new Polygonizer(true);
399 Geometry fixed
= GeometryFixer
.fix(polygon
, false);
400 polygonizer
.add(fixed
);
401 @SuppressWarnings("unchecked")
402 List
<Polygon
> polygons
= new ArrayList
<>(polygonizer
.getPolygons());
404 if (polygons
.size() == 0) {
405 throw new IllegalStateException("Polygonizer failed to extract any polygon");
409 if (polygons
.size() == 1) {
410 best
= polygons
.get(0);
412 double totalArea
= fixed
.getArea();
413 best
= polygons
.get(0);
414 double bestArea
= best
.getArea();
415 for (int i
= 1; i
< polygons
.size(); i
++) {
416 Polygon p
= polygons
.get(i
);
417 double a
= p
.getArea();
422 // double ratio = a / totalArea;
425 double ignoredRatio
= (totalArea
- bestArea
) / totalArea
;
426 if (ignoredRatio
> ignoredAreaRatio
)
429 if (!best
.isValid() || !best
.isSimple()) {
430 throw new IllegalStateException("The polygon is not simple and/or valid after cleaning");
433 // while we are here, we make sure that the geometry will be normalised
439 * The smallest polygon without holes containing all the points in this
442 public static Polygon
concaveHull(Geometry geom
, double lengthRatio
) {
443 Objects
.requireNonNull(geom
);
444 if (geom
.getNumPoints() < 3)
445 throw new IllegalArgumentException("At least 3 points are reuired to compute the concave hull and geometry "
446 + geom
.getClass() + " contains only " + geom
.getNumPoints());
447 Geometry hull
= ConcaveHull
.concaveHullByLengthRatio(geom
, lengthRatio
, false);
448 if (hull
instanceof Polygon polygon
)
451 throw new IllegalStateException("Hull is not a polygon but a " + hull
.getClass());