org.argeo.app.geo/src/org/argeo/app/geo/GeoUtils.java

   1 package org.argeo.app.geo;
   2
   3 import java.io.IOException;
   4 import java.io.Serializable;
   5 import java.io.Writer;
   6 import java.nio.file.Path;
   7 import java.util.ArrayList;
   8 import java.util.HashMap;
   9 import java.util.List;
  10 import java.util.Map;
  11 import java.util.Objects;
  12
  13 import javax.measure.Quantity;
  14 import javax.measure.quantity.Area;
  15 import javax.xml.transform.TransformerException;
  16
  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;
  65
  66 import si.uom.SI;
  67 import tech.units.indriya.quantity.Quantities;
  68
  69 /** Utilities around geographical format, mostly wrapping GeoTools patterns. */
  70 public class GeoUtils {
  71         public final static String EPSG_4326 = "EPSG:4326";
  72
  73         /** In square meters. */
  74         public static Quantity<Area> calcArea(SimpleFeature feature) {
  75                 try {
  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");
  83
  84                         // create automatic CRS for optimal computation
  85                         String code = "AUTO:42001," + centroid.getX() + "," + centroid.getY();
  86                         CoordinateReferenceSystem auto = CRS.decode(code);
  87
  88                         MathTransform transform = CRS.findMathTransform(DefaultGeographicCRS.WGS84, auto);
  89
  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);
  94                 }
  95         }
  96
  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());
 101         }
 102
 103         public static Quantity<Area> calcArea(Polygon polygon, CoordinateReferenceSystem crs, Point wgs84centroid) {
 104                 try {
 105                         // create automatic CRS for optimal computation
 106                         String code = "AUTO:42001," + wgs84centroid.getX() + "," + wgs84centroid.getY();
 107                         CoordinateReferenceSystem auto = CRS.decode(code);
 108
 109                         MathTransform transform = CRS.findMathTransform(crs, auto);
 110
 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);
 115                 }
 116         }
 117
 118         public static void exportToSvg(Geometry[] geometries, Writer out, int width, int height) {
 119                 try {
 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=\"");
 128
 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);
 133
 134                                         MathTransform transform = CRS.findMathTransform(DefaultGeographicCRS.WGS84, auto);
 135
 136                                         Polygon projed = (Polygon) JTS.transform(p, transform);
 137
 138                                         // EPSG4326 are in lat/lon order, so we invert coordinates
 139                                         for (Coordinate coord : projed.getCoordinates()) {
 140                                                 double x = coord.y;
 141                                                 if (x < minX)
 142                                                         minX = x;
 143                                                 if (x > maxX)
 144                                                         maxX = x;
 145                                                 double y = -coord.x;
 146                                                 if (y < minY)
 147                                                         minY = y;
 148                                                 if (y > maxY)
 149                                                         maxY = y;
 150                                                 sb.append(x + "," + y + " ");
 151                                         }
 152                                         sb.append("\">");
 153                                         sb.append("</polyline>\n");
 154                                         shapes.add(sb.toString());
 155                                 } else {
 156                                         throw new IllegalArgumentException("Unsuppported geometry type " + geometry.getClass().getName());
 157                                 }
 158                         }
 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");
 166                         out.write(">\n");
 167                         for (String shape : shapes) {
 168                                 out.write(shape);
 169                                 out.write("\n");
 170                         }
 171                         out.write("</svg>");
 172                 } catch (IOException | FactoryException | MismatchedDimensionException | TransformException e) {
 173                         throw new RuntimeException("Cannot export to SVG", e);
 174                 }
 175         }
 176
 177         /** Write a list of simple features to a shapefile. */
 178         public static void saveFeaturesAsShapefile(SimpleFeatureType featureType, List<SimpleFeature> features,
 179                         Path shpFile) {
 180                 try {
 181                         ShapefileDataStoreFactory dataStoreFactory = new ShapefileDataStoreFactory();
 182
 183                         Map<String, Serializable> params = new HashMap<>();
 184                         params.put("url", shpFile.toUri().toURL());
 185
 186                         params.put("create spatial index", Boolean.TRUE);
 187
 188                         ShapefileDataStore newDataStore = (ShapefileDataStore) dataStoreFactory.createNewDataStore(params);
 189                         newDataStore.createSchema(featureType);
 190
 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);
 196
 197                                 try (Transaction transaction = new DefaultTransaction("create")) {
 198                                         try {
 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);
 205                                         }
 206                                 }
 207                         } else {
 208                                 throw new IllegalArgumentException(typeName + " does not support read/write access");
 209                         }
 210                 } catch (IOException e) {
 211                         throw new RuntimeException("Cannot write shapefile " + shpFile, e);
 212                 }
 213         }
 214
 215         /*
 216          * STYLING
 217          */
 218
 219         public static Style createStyleFromCss(String css) {
 220                 Stylesheet ss = CssParser.parse(css);
 221                 CssTranslator translator = new CssTranslator();
 222                 Style style = translator.translate(ss);
 223
 224 //              try {
 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();
 231 //              }
 232
 233                 return (Style) style;
 234         }
 235
 236         public static String createSldFromCss(String name, String title, String css) {
 237                 StyledLayerDescriptor sld = GeoTools.STYLE_FACTORY.createStyledLayerDescriptor();
 238                 sld.setName(name);
 239                 sld.setTitle(title);
 240
 241                 NamedLayer layer = GeoTools.STYLE_FACTORY.createNamedLayer();
 242                 layer.setName(name);
 243
 244                 Style style = createStyleFromCss(css);
 245                 layer.styles().add(style);
 246
 247                 sld.layers().add(layer);
 248                 return sldToXml(sld);
 249         }
 250
 251         public static String sldToXml(StyledLayerDescriptor sld) {
 252                 try {
 253                         SLDTransformer styleTransform = new SLDTransformer();
 254                         String xml = styleTransform.transform(sld);
 255 //                      System.out.println(xml);
 256                         return xml;
 257                 } catch (TransformerException e) {
 258                         throw new IllegalStateException("Cannot write SLD as XML", e);
 259                 }
 260         }
 261
 262         public static void main(String... args) {
 263                 String css = """
 264                                 * {
 265                                    mark: symbol(circle);
 266                                    mark-size: 6px;
 267                                  }
 268
 269                                  :mark {
 270                                    fill: red;
 271                                  }
 272
 273                                                                 """;
 274                 createSldFromCss("test", "Test", css);
 275         }
 276
 277         public static String createTestSLD() {
 278
 279                 StyleFactory sf = CommonFactoryFinder.getStyleFactory();
 280                 FilterFactory ff = CommonFactoryFinder.getFilterFactory();
 281
 282                 StyledLayerDescriptor sld = sf.createStyledLayerDescriptor();
 283                 sld.setName("sld");
 284                 sld.setTitle("Example");
 285                 sld.setAbstract("Example Style Layer Descriptor");
 286
 287                 UserLayer layer = sf.createUserLayer();
 288                 layer.setName("layer");
 289
 290                 //
 291                 // define constraint limited what features the sld applies to
 292                 FeatureTypeConstraint constraint = sf.createFeatureTypeConstraint("Feature", Filter.INCLUDE);
 293
 294                 layer.layerFeatureConstraints().add(constraint);
 295
 296                 //
 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");
 302
 303                 //
 304                 // define feature type styles used to actually define how features are rendered
 305                 FeatureTypeStyle featureTypeStyle = sf.createFeatureTypeStyle();
 306
 307                 // RULE 1
 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)));
 314
 315                 //
 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));
 319
 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
 323 //
 324 //        OnLineResourceImpl png = new OnLineResourceImpl(new URI("file:city.png"));
 325 //        png.freeze(); // freeze to prevent modification at runtime
 326
 327                 //
 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
 334
 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
 340
 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);
 344
 345                 rule1.symbolizers().add(pointSymbolizer);
 346
 347                 featureTypeStyle.rules().add(rule1);
 348
 349                 //
 350                 // RULE 2 Default
 351
 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);
 361
 362                 style.featureTypeStyles().add(featureTypeStyle);
 363
 364                 layer.userStyles().add(style);
 365
 366                 sld.layers().add(layer);
 367
 368                 try {
 369                         SLDTransformer styleTransform = new SLDTransformer();
 370                         String xml = styleTransform.transform(sld);
 371                         System.out.println(xml);
 372                         return xml;
 373                 } catch (TransformerException e) {
 374                         throw new IllegalStateException(e);
 375                 }
 376
 377         }
 378
 379         public static long getScaleFromResolution(long resolution) {
 380                 // see
 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);
 386         }
 387
 388         /*
 389          * GEOMETRY
 390          */
 391         /**
 392          * Ensure that a {@link Polygon} is valid and simple by removing artefacts
 393          * (typically coming from GPS).
 394          *
 395          * @return a simple and valid polygon, or null if the ignoredArea ratio is above
 396          *         the provided threshold.
 397          */
 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());
 404
 405                 if (polygons.size() == 0) {
 406                         throw new IllegalStateException("Polygonizer failed to extract any polygon");
 407                 }
 408
 409                 Polygon best;
 410                 if (polygons.size() == 1) {
 411                         best = polygons.get(0);
 412                 } else {
 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();
 419                                 if (a > bestArea) {
 420                                         best = p;
 421                                         bestArea = a;
 422                                 } else {
 423                                         // double ratio = a / totalArea;
 424                                 }
 425                         }
 426                         double ignoredRatio = (totalArea - bestArea) / totalArea;
 427                         if (ignoredRatio > ignoredAreaRatio)
 428                                 return null;
 429
 430                         if (!best.isValid() || !best.isSimple()) {
 431                                 throw new IllegalStateException("The polygon is not simple and/or valid after cleaning");
 432                         }
 433                 }
 434                 // while we are here, we make sure that the geometry will be normalised
 435                 best.normalize();
 436                 return best;
 437         }
 438
 439         /**
 440          * The smallest polygon without holes containing all the points in this
 441          * geometry.
 442          */
 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)
 450                         return polygon;
 451                 else
 452                         throw new IllegalStateException("Hull is not a polygon but a " + hull.getClass());
 453         }
 454
 455         /** Singleton. */
 456         private GeoUtils() {
 457         }
 458 }