package com.onthegomap.planetiler.geo; import com.onthegomap.planetiler.collection.LongLongMap; import com.onthegomap.planetiler.config.PlanetilerConfig; import com.onthegomap.planetiler.stats.Stats; import java.util.ArrayList; import java.util.List; import java.util.stream.Stream; import org.locationtech.jts.algorithm.Area; import org.locationtech.jts.geom.Coordinate; import org.locationtech.jts.geom.CoordinateSequence; import org.locationtech.jts.geom.CoordinateXY; import org.locationtech.jts.geom.Envelope; import org.locationtech.jts.geom.Geometry; import org.locationtech.jts.geom.GeometryCollection; import org.locationtech.jts.geom.GeometryFactory; import org.locationtech.jts.geom.LineSegment; import org.locationtech.jts.geom.LineString; import org.locationtech.jts.geom.LinearRing; import org.locationtech.jts.geom.MultiPolygon; import org.locationtech.jts.geom.Point; import org.locationtech.jts.geom.Polygon; import org.locationtech.jts.geom.PrecisionModel; import org.locationtech.jts.geom.TopologyException; import org.locationtech.jts.geom.impl.PackedCoordinateSequence; import org.locationtech.jts.geom.impl.PackedCoordinateSequenceFactory; import org.locationtech.jts.geom.util.GeometryFixer; import org.locationtech.jts.geom.util.GeometryTransformer; import org.locationtech.jts.io.WKBReader; import org.locationtech.jts.precision.GeometryPrecisionReducer; /** * A collection of utilities for working with JTS data structures and geographic data. *
* "world" coordinates in this class refer to web mercator coordinates where the top-left/northwest corner of the map is
* (0,0) and bottom-right/southeast corner is (1,1).
*/
public class GeoUtils {
/** Rounding precision for 256x256px tiles encoded using 4096 values. */
public static final PrecisionModel TILE_PRECISION = new PrecisionModel(4096d / 256d);
public static final GeometryFactory JTS_FACTORY = new GeometryFactory(PackedCoordinateSequenceFactory.DOUBLE_FACTORY);
public static final WKBReader WKB_READER = new WKBReader(JTS_FACTORY);
public static final Geometry EMPTY_GEOMETRY = JTS_FACTORY.createGeometryCollection();
public static final Point EMPTY_POINT = JTS_FACTORY.createPoint();
public static final LineString EMPTY_LINE = JTS_FACTORY.createLineString();
public static final Polygon EMPTY_POLYGON = JTS_FACTORY.createPolygon();
private static final LineString[] EMPTY_LINE_STRING_ARRAY = new LineString[0];
private static final Polygon[] EMPTY_POLYGON_ARRAY = new Polygon[0];
private static final Point[] EMPTY_POINT_ARRAY = new Point[0];
private static final double WORLD_RADIUS_METERS = 6_378_137;
public static final double WORLD_CIRCUMFERENCE_METERS = Math.PI * 2 * WORLD_RADIUS_METERS;
private static final double RADIANS_PER_DEGREE = Math.PI / 180;
private static final double DEGREES_PER_RADIAN = 180 / Math.PI;
private static final double LOG2 = Math.log(2);
/**
* Transform web mercator coordinates where top-left corner of the planet is (0,0) and bottom-right is (1,1) to
* latitude/longitude coordinates.
*/
private static final GeometryTransformer UNPROJECT_WORLD_COORDS = new GeometryTransformer() {
@Override
protected CoordinateSequence transformCoordinates(CoordinateSequence coords, Geometry parent) {
CoordinateSequence copy = new PackedCoordinateSequence.Double(coords.size(), 2, 0);
for (int i = 0; i < coords.size(); i++) {
copy.setOrdinate(i, 0, getWorldLon(coords.getX(i)));
copy.setOrdinate(i, 1, getWorldLat(coords.getY(i)));
}
return copy;
}
};
/**
* Transform latitude/longitude coordinates to web mercator where top-left corner of the planet is (0,0) and
* bottom-right is (1,1).
*/
private static final GeometryTransformer PROJECT_WORLD_COORDS = new GeometryTransformer() {
@Override
protected CoordinateSequence transformCoordinates(CoordinateSequence coords, Geometry parent) {
CoordinateSequence copy = new PackedCoordinateSequence.Double(coords.size(), 2, 0);
for (int i = 0; i < coords.size(); i++) {
copy.setOrdinate(i, 0, getWorldX(coords.getX(i)));
copy.setOrdinate(i, 1, getWorldY(coords.getY(i)));
}
return copy;
}
};
private static final double MAX_LAT = getWorldLat(-0.1);
private static final double MIN_LAT = getWorldLat(1.1);
// to pack latitude/longitude into a single long, we round them to 31 bits of precision
private static final double QUANTIZED_WORLD_SIZE = Math.pow(2, 31);
private static final double HALF_QUANTIZED_WORLD_SIZE = QUANTIZED_WORLD_SIZE / 2;
private static final long LOWER_32_BIT_MASK = (1L << 32) - 1L;
/**
* Bounds for the entire area of the planet that a web mercator projection covers, where top left is (0,0) and bottom
* right is (1,1).
*/
public static final Envelope WORLD_BOUNDS = new Envelope(0, 1, 0, 1);
/**
* Bounds for the entire area of the planet that a web mercator projection covers in latitude/longitude coordinates.
*/
public static final Envelope WORLD_LAT_LON_BOUNDS = toLatLonBoundsBounds(WORLD_BOUNDS);
// should not instantiate
private GeoUtils() {}
/**
* Returns a copy of {@code geom} transformed from latitude/longitude coordinates to web mercator where top-left
* corner of the planet is (0,0) and bottom-right is (1,1).
*/
public static Geometry latLonToWorldCoords(Geometry geom) {
return PROJECT_WORLD_COORDS.transform(geom);
}
/**
* Returns a copy of {@code geom} transformed from web mercator where top-left corner of the planet is (0,0) and
* bottom-right is (1,1) to latitude/longitude.
*/
public static Geometry worldToLatLonCoords(Geometry geom) {
return UNPROJECT_WORLD_COORDS.transform(geom);
}
/**
* Returns a copy of {@code worldBounds} transformed from web mercator where top-left corner of the planet is (0,0)
* and bottom-right is (1,1) to latitude/longitude.
*/
public static Envelope toLatLonBoundsBounds(Envelope worldBounds) {
return new Envelope(
getWorldLon(worldBounds.getMinX()),
getWorldLon(worldBounds.getMaxX()),
getWorldLat(worldBounds.getMinY()),
getWorldLat(worldBounds.getMaxY()));
}
/**
* Returns a copy of {@code lonLatBounds} transformed from latitude/longitude coordinates to web mercator where
* top-left corner of the planet is (0,0) and bottom-right is (1,1).
*/
public static Envelope toWorldBounds(Envelope lonLatBounds) {
return new Envelope(
getWorldX(lonLatBounds.getMinX()),
getWorldX(lonLatBounds.getMaxX()),
getWorldY(lonLatBounds.getMinY()),
getWorldY(lonLatBounds.getMaxY())
);
}
/**
* Returns the longitude for a web mercator coordinate {@code x} where 0 is the international date line on the west
* side, 1 is the international date line on the east side, and 0.5 is the prime meridian.
*/
public static double getWorldLon(double x) {
return x * 360 - 180;
}
/**
* Returns the latitude for a web mercator {@code y} coordinate where 0 is the north edge of the map, 0.5 is the
* equator, and 1 is the south edge of the map.
*/
public static double getWorldLat(double y) {
double n = Math.PI - 2 * Math.PI * y;
return DEGREES_PER_RADIAN * Math.atan(0.5 * (Math.exp(n) - Math.exp(-n)));
}
/**
* Returns the web mercator X coordinate for {@code longitude} where 0 is the international date line on the west
* side, 1 is the international date line on the east side, and 0.5 is the prime meridian.
*/
public static double getWorldX(double longitude) {
return (longitude + 180) / 360;
}
/**
* Returns the web mercator Y coordinate for {@code latitude} where 0 is the north edge of the map, 0.5 is the
* equator, and 1 is the south edge of the map.
*/
public static double getWorldY(double latitude) {
if (latitude <= MIN_LAT) {
return 1.1;
}
if (latitude >= MAX_LAT) {
return -0.1;
}
double sin = Math.sin(latitude * RADIANS_PER_DEGREE);
return 0.5 - 0.25 * Math.log((1 + sin) / (1 - sin)) / Math.PI;
}
/**
* Returns a latitude/longitude coordinate encoded into a single 64-bit long for storage in a {@link LongLongMap} that
* can be decoded using {@link #decodeWorldX(long)} and {@link #decodeWorldY(long)}.
*/
public static long encodeFlatLocation(double lon, double lat) {
double worldX = getWorldX(lon) + 1;
double worldY = getWorldY(lat) + 1;
long x = (long) (worldX * HALF_QUANTIZED_WORLD_SIZE);
long y = (long) (worldY * HALF_QUANTIZED_WORLD_SIZE);
return (x << 32) | (y & LOWER_32_BIT_MASK);
}
/**
* Returns the web mercator Y coordinate of the latitude/longitude encoded with
* {@link #encodeFlatLocation(double, double)}.
*/
public static double decodeWorldY(long encoded) {
return ((encoded & LOWER_32_BIT_MASK) / HALF_QUANTIZED_WORLD_SIZE) - 1;
}
/**
* Returns the web mercator X coordinate of the latitude/longitude encoded with
* {@link #encodeFlatLocation(double, double)}.
*/
public static double decodeWorldX(long encoded) {
return ((encoded >>> 32) / HALF_QUANTIZED_WORLD_SIZE) - 1;
}
/**
* Returns an approximate zoom level that a map should be displayed at to show all of {@code envelope}, specified in
* latitude/longitude coordinates.
*/
public static double getZoomFromLonLatBounds(Envelope envelope) {
Envelope worldBounds = GeoUtils.toWorldBounds(envelope);
return getZoomFromWorldBounds(worldBounds);
}
/**
* Returns an approximate zoom level that a map should be displayed at to show all of {@code envelope}, specified in
* web mercator coordinates.
*/
public static double getZoomFromWorldBounds(Envelope worldBounds) {
double maxEdge = Math.max(worldBounds.getWidth(), worldBounds.getHeight());
return Math.max(0, -Math.log(maxEdge) / Math.log(2));
}
/** Returns the width in meters of a single pixel of a 256x256 px tile at the given {@code zoom} level. */
public static double metersPerPixelAtEquator(int zoom) {
return WORLD_CIRCUMFERENCE_METERS / Math.pow(2d, zoom + 8d);
}
/** Returns the length in pixels for a given number of meters on a 256x256 px tile at the given {@code zoom} level. */
public static double metersToPixelAtEquator(int zoom, double meters) {
return meters / metersPerPixelAtEquator(zoom);
}
public static Point point(double x, double y) {
return JTS_FACTORY.createPoint(new CoordinateXY(x, y));
}
public static Point point(Coordinate coord) {
return JTS_FACTORY.createPoint(coord);
}
public static Geometry createMultiLineString(List
* The result will be clamped to the range [0, {@link PlanetilerConfig#MAX_MAXZOOM}].
*/
public static int minZoomForPixelSize(double worldGeometrySize, double minPixelSize) {
double worldPixels = worldGeometrySize * 256;
return Math.clamp((int) Math.ceil(Math.log(minPixelSize / worldPixels) / LOG2), 0,
PlanetilerConfig.MAX_MAXZOOM);
}
/** Helper class to sort polygons by area of their outer shell. */
private record PolyAndArea(Polygon poly, double area) implements Comparable