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                                         for (Coordinate coord : projed.getCoordinates()) {
 139                                                 double x = coord.x;
 140                                                 if (x < minX)
 141                                                         minX = x;
 142                                                 if (x > maxX)
 143                                                         maxX = x;
 144                                                 double y = -coord.y;
 145                                                 if (y < minY)
 146                                                         minY = y;
 147                                                 if (y > maxY)
 148                                                         maxY = y;
 149                                                 sb.append(x + "," + y + " ");
 150                                         }
 151                                         sb.append("\">");
 152                                         sb.append("</polyline>\n");
 153                                         shapes.add(sb.toString());
 154                                 } else {
 155                                         throw new IllegalArgumentException("Unsuppported geometry type " + geometry.getClass().getName());
 156                                 }
 157                         }
 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");
 165                         out.write(">\n");
 166                         for (String shape : shapes) {
 167                                 out.write(shape);
 168                                 out.write("\n");
 169                         }
 170                         out.write("</svg>");
 171                 } catch (IOException | FactoryException | MismatchedDimensionException | TransformException e) {
 172                         throw new RuntimeException("Cannot export to SVG", e);
 173                 }
 174         }
 175
 176         /** Write a list of simple features to a shapefile. */
 177         public static void saveFeaturesAsShapefile(SimpleFeatureType featureType, List<SimpleFeature> features,
 178                         Path shpFile) {
 179                 try {
 180                         ShapefileDataStoreFactory dataStoreFactory = new ShapefileDataStoreFactory();
 181
 182                         Map<String, Serializable> params = new HashMap<>();
 183                         params.put("url", shpFile.toUri().toURL());
 184
 185                         params.put("create spatial index", Boolean.TRUE);
 186
 187                         ShapefileDataStore newDataStore = (ShapefileDataStore) dataStoreFactory.createNewDataStore(params);
 188                         newDataStore.createSchema(featureType);
 189
 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);
 195
 196                                 try (Transaction transaction = new DefaultTransaction("create")) {
 197                                         try {
 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);
 204                                         }
 205                                 }
 206                         } else {
 207                                 throw new IllegalArgumentException(typeName + " does not support read/write access");
 208                         }
 209                 } catch (IOException e) {
 210                         throw new RuntimeException("Cannot write shapefile " + shpFile, e);
 211                 }
 212         }
 213
 214         /*
 215          * STYLING
 216          */
 217
 218         public static Style createStyleFromCss(String css) {
 219                 Stylesheet ss = CssParser.parse(css);
 220                 CssTranslator translator = new CssTranslator();
 221                 Style style = translator.translate(ss);
 222
 223 //              try {
 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();
 230 //              }
 231
 232                 return (Style) style;
 233         }
 234
 235         public static String createSldFromCss(String name, String title, String css) {
 236                 StyledLayerDescriptor sld = GeoTools.STYLE_FACTORY.createStyledLayerDescriptor();
 237                 sld.setName(name);
 238                 sld.setTitle(title);
 239
 240                 NamedLayer layer = GeoTools.STYLE_FACTORY.createNamedLayer();
 241                 layer.setName(name);
 242
 243                 Style style = createStyleFromCss(css);
 244                 layer.styles().add(style);
 245
 246                 sld.layers().add(layer);
 247                 return sldToXml(sld);
 248         }
 249
 250         public static String sldToXml(StyledLayerDescriptor sld) {
 251                 try {
 252                         SLDTransformer styleTransform = new SLDTransformer();
 253                         String xml = styleTransform.transform(sld);
 254 //                      System.out.println(xml);
 255                         return xml;
 256                 } catch (TransformerException e) {
 257                         throw new IllegalStateException("Cannot write SLD as XML", e);
 258                 }
 259         }
 260
 261         public static void main(String... args) {
 262                 String css = """
 263                                 * {
 264                                    mark: symbol(circle);
 265                                    mark-size: 6px;
 266                                  }
 267
 268                                  :mark {
 269                                    fill: red;
 270                                  }
 271
 272                                                                 """;
 273                 createSldFromCss("test", "Test", css);
 274         }
 275
 276         public static String createTestSLD() {
 277
 278                 StyleFactory sf = CommonFactoryFinder.getStyleFactory();
 279                 FilterFactory ff = CommonFactoryFinder.getFilterFactory();
 280
 281                 StyledLayerDescriptor sld = sf.createStyledLayerDescriptor();
 282                 sld.setName("sld");
 283                 sld.setTitle("Example");
 284                 sld.setAbstract("Example Style Layer Descriptor");
 285
 286                 UserLayer layer = sf.createUserLayer();
 287                 layer.setName("layer");
 288
 289                 //
 290                 // define constraint limited what features the sld applies to
 291                 FeatureTypeConstraint constraint = sf.createFeatureTypeConstraint("Feature", Filter.INCLUDE);
 292
 293                 layer.layerFeatureConstraints().add(constraint);
 294
 295                 //
 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");
 301
 302                 //
 303                 // define feature type styles used to actually define how features are rendered
 304                 FeatureTypeStyle featureTypeStyle = sf.createFeatureTypeStyle();
 305
 306                 // RULE 1
 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)));
 313
 314                 //
 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));
 318
 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
 322 //
 323 //        OnLineResourceImpl png = new OnLineResourceImpl(new URI("file:city.png"));
 324 //        png.freeze(); // freeze to prevent modification at runtime
 325
 326                 //
 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
 333
 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
 339
 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);
 343
 344                 rule1.symbolizers().add(pointSymbolizer);
 345
 346                 featureTypeStyle.rules().add(rule1);
 347
 348                 //
 349                 // RULE 2 Default
 350
 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);
 360
 361                 style.featureTypeStyles().add(featureTypeStyle);
 362
 363                 layer.userStyles().add(style);
 364
 365                 sld.layers().add(layer);
 366
 367                 try {
 368                         SLDTransformer styleTransform = new SLDTransformer();
 369                         String xml = styleTransform.transform(sld);
 370                         System.out.println(xml);
 371                         return xml;
 372                 } catch (TransformerException e) {
 373                         throw new IllegalStateException(e);
 374                 }
 375
 376         }
 377
 378         public static long getScaleFromResolution(long resolution) {
 379                 // see
 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);
 385         }
 386
 387         /*
 388          * GEOMETRY
 389          */
 390         /**
 391          * Ensure that a {@link Polygon} is valid and simple by removing artefacts
 392          * (typically coming from GPS).
 393          *
 394          * @return a simple and valid polygon, or null if the ignoredArea ratio is above
 395          *         the provided threshold.
 396          */
 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());
 403
 404                 if (polygons.size() == 0) {
 405                         throw new IllegalStateException("Polygonizer failed to extract any polygon");
 406                 }
 407
 408                 Polygon best;
 409                 if (polygons.size() == 1) {
 410                         best = polygons.get(0);
 411                 } else {
 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();
 418                                 if (a > bestArea) {
 419                                         best = p;
 420                                         bestArea = a;
 421                                 } else {
 422                                         // double ratio = a / totalArea;
 423                                 }
 424                         }
 425                         double ignoredRatio = (totalArea - bestArea) / totalArea;
 426                         if (ignoredRatio > ignoredAreaRatio)
 427                                 return null;
 428
 429                         if (!best.isValid() || !best.isSimple()) {
 430                                 throw new IllegalStateException("The polygon is not simple and/or valid after cleaning");
 431                         }
 432                 }
 433                 // while we are here, we make sure that the geometry will be normalised
 434                 best.normalize();
 435                 return best;
 436         }
 437
 438         /**
 439          * The smallest polygon without holes containing all the points in this
 440          * geometry.
 441          */
 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)
 449                         return polygon;
 450                 else
 451                         throw new IllegalStateException("Hull is not a polygon but a " + hull.getClass());
 452         }
 453
 454         /** Singleton. */
 455         private GeoUtils() {
 456         }
 457 }