第19章 - 坐标系、单位与工程陷阱

import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.geom.GeometryFactory;
import org.locationtech.jts.geom.PrecisionModel;
import org.locationtech.jts.io.WKTReader;

public class WktSridDemo {
    public static void main(String[] args) throws Exception {
        GeometryFactory gf4326 =
            new GeometryFactory(new PrecisionModel(PrecisionModel.FLOATING), 4326);
        WKTReader reader = new WKTReader(gf4326);

        Geometry g = reader.read("POINT (116.4074 39.9042)");
        System.out.println(g.getSRID());   // 4326（来自工厂，不是 WKT）

        WKTReader defaultReader = new WKTReader(); // 默认工厂 SRID = 0
        Geometry g0 = defaultReader.read("POINT (116.4074 39.9042)");
        System.out.println(g0.getSRID());  // 0

        // 解析 EWKT 前缀会失败（JTS 标准 reader 不支持）
        try {
            defaultReader.read("SRID=4326;POINT (116.4074 39.9042)");
        } catch (Exception ex) {
            System.out.println("EWKT 前缀不被识别: " + ex.getClass().getSimpleName());
        }
    }
}

import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.geom.GeometryFactory;
import org.locationtech.jts.geom.PrecisionModel;
import org.locationtech.jts.io.WKBReader;
import org.locationtech.jts.io.WKBWriter;

public class WkbSridDemo {
    public static void main(String[] args) throws Exception {
        GeometryFactory gf =
            new GeometryFactory(new PrecisionModel(), 4490);
        Geometry g = gf.createPoint(new org.locationtech.jts.geom.Coordinate(116.4, 39.9));

        // 1) 标准 WKB：SRID 丢失
        byte[] plain = new WKBWriter().write(g);
        Geometry back1 = new WKBReader().read(plain);
        System.out.println(back1.getSRID());       // 0（默认工厂）

        // 2) EWKB：SRID 保留
        byte[] ewkb = new WKBWriter(2, true).write(g);
        Geometry back2 = new WKBReader().read(ewkb);
        System.out.println(back2.getSRID());       // 4490
    }
}

import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.io.geojson.GeoJsonReader;
import org.locationtech.jts.io.geojson.GeoJsonWriter;

public class GeoJsonSridDemo {
    public static void main(String[] args) throws Exception {
        GeoJsonReader reader = new GeoJsonReader();
        Geometry g = reader.read(
            "{\"type\":\"Point\",\"coordinates\":[116.4074,39.9042]}");
        System.out.println(g.getSRID());   // 4326（GeoJSON 默认）

        // 若数据实际是投影坐标，读取后必须显式改标签（但坐标不会转换！）
        g.setSRID(4547);
        System.out.println(g.getSRID());   // 4547

        GeoJsonWriter writer = new GeoJsonWriter();
        System.out.println(writer.write(g).contains("crs"));
        // 是否输出 crs 取决于 SRID 与 writer 配置
    }
}

import org.geotools.geometry.jts.JTS;
import org.geotools.referencing.CRS;
import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.io.WKTReader;
import org.opengis.referencing.crs.CoordinateReferenceSystem;
import org.opengis.referencing.operation.MathTransform;

public class GeoToolsTransformDemo {
    public static void main(String[] args) throws Exception {
        WKTReader reader = new WKTReader();
        Geometry lonLat = reader.read("POINT (116.4074 39.9042)");

        // true = 强制经度在前，避免 lat/lon 轴顺序问题
        CoordinateReferenceSystem src = CRS.decode("EPSG:4326", true);
        CoordinateReferenceSystem dst = CRS.decode("EPSG:4547", true);

        MathTransform mt = CRS.findMathTransform(src, dst, true);
        Geometry projected = JTS.transform(lonLat, mt);
        System.out.println(projected);            // 投影后的米制坐标
        System.out.println(projected.getSRID());  // 0：JTS.transform 不会自动设置 SRID

        // 逆变换回经纬度
        MathTransform inverse = mt.inverse();
        Geometry back = JTS.transform(projected, inverse);
        System.out.println(back);
    }
}

import org.locationtech.proj4j.CRSFactory;
import org.locationtech.proj4j.CoordinateReferenceSystem;
import org.locationtech.proj4j.CoordinateTransform;
import org.locationtech.proj4j.CoordinateTransformFactory;
import org.locationtech.proj4j.ProjCoordinate;

public class Proj4jDemo {
    public static void main(String[] args) {
        CRSFactory crsFactory = new CRSFactory();
        CoordinateReferenceSystem src = crsFactory.createFromName("EPSG:4326");
        CoordinateReferenceSystem dst = crsFactory.createFromName("EPSG:3857");

        CoordinateTransformFactory ctFactory = new CoordinateTransformFactory();
        CoordinateTransform transform = ctFactory.createTransform(src, dst);

        ProjCoordinate in = new ProjCoordinate(116.4074, 39.9042);
        ProjCoordinate out = new ProjCoordinate();
        transform.transform(in, out);
        System.out.println(out.x + ", " + out.y);
        // 约 1.296e7, 4.85e6（Web Mercator 米）
    }
}

import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.io.WKTReader;

public class DegreeAreaDemo {
    public static void main(String[] args) throws Exception {
        WKTReader r = new WKTReader();
        // 北纬 40~41 度、东经 116~117 度的一个方块
        Geometry box = r.read(
            "POLYGON ((116 40, 117 40, 117 41, 116 41, 116 40))");

        double areaDeg2 = box.getArea();
        System.out.println("平面面积 = " + areaDeg2 + " 平方度"); // 1.0

        double latMid = Math.toRadians(40.5);
        double km2PerDeg2 = 111.32 * 110.57 * Math.cos(latMid);
        System.out.println("近似面积 = " + (areaDeg2 * km2PerDeg2) + " 平方公里");
        // 约 9365 平方公里（球面精确值约 9401）
    }
}

import org.locationtech.jts.geom.Coordinate;
import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.geom.GeometryFactory;

public class HaversineDemo {
    static final double R = 6371008.8; // 平均地球半径，米

    public static double haversine(Coordinate p1, Coordinate p2) {
        double lat1 = Math.toRadians(p1.y);
        double lat2 = Math.toRadians(p2.y);
        double dLat = lat2 - lat1;
        double dLon = Math.toRadians(p2.x - p1.x);
        double a = Math.sin(dLat / 2) * Math.sin(dLat / 2)
                 + Math.cos(lat1) * Math.cos(lat2)
                 * Math.sin(dLon / 2) * Math.sin(dLon / 2);
        double c = 2 * Math.asin(Math.sqrt(a));
        return R * c;
    }

    public static void main(String[] args) {
        Coordinate beijing = new Coordinate(116.4074, 39.9042);
        Coordinate shanghai = new Coordinate(121.4737, 31.2304);

        GeometryFactory gf = new GeometryFactory();
        Geometry line = gf.createLineString(new Coordinate[] {beijing, shanghai});

        System.out.println("度平面长度 = " + line.getLength() + " 度");
        System.out.println("Haversine = " + haversine(beijing, shanghai) / 1000 + " km");
        System.out.println("度乘 111.32 = " + line.getLength() * 111.32 + " km");
    }
}

import org.locationtech.jts.geom.Coordinate;
import org.locationtech.jts.geom.Envelope;
import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.geom.GeometryFactory;
import org.locationtech.jts.geom.util.AffineTransformation;
import org.locationtech.jts.io.WKTReader;

public class AntimeridianDemo {
    // 1) 逐坐标归一化：负经度加 360，统一到 [0, 360)
    static Geometry normalize0360(Geometry g) {
        Geometry copy = g.copy();
        copy.apply(new org.locationtech.jts.geom.CoordinateFilter() {
            public void filter(Coordinate c) {
                if (c.x < 0) c.x += 360;
            }
        });
        return copy;
    }

    public static Geometry[] split(Geometry crossing) {
        GeometryFactory gf = crossing.getFactory();
        Geometry shifted = normalize0360(crossing);

        // 2) 用两个半球矩形切分
        Geometry west = shifted.intersection(
            gf.toGeometry(new Envelope(0, 180, -90, 90)));   // 0~180
        Geometry east = shifted.intersection(
            gf.toGeometry(new Envelope(180, 360, -90, 90))); // 180~360

        // 3) 把东半部分移回 [-180, 0]
        if (!east.isEmpty()) {
            east = AffineTransformation.translationInstance(-360, 0).transform(east);
        }
        return new Geometry[] {west, east};
    }

    public static void main(String[] args) throws Exception {
        WKTReader r = new WKTReader();
        Geometry island = r.read(
            "POLYGON ((179 10, -179 10, -179 12, 179 12, 179 10))");

        System.out.println(island.getEnvelopeInternal().getWidth()); // 358.0

        Geometry[] parts = split(island);
        System.out.println(parts[0]);
        System.out.println(parts[1]);
        System.out.println(parts[0].getArea() + parts[1].getArea()); // 4.0 平方度（2×2）
    }
}

public class MercatorDistortionDemo {
    static final double R = 6378137.0;

    public static double[] forward(double lonDeg, double latDeg) {
        double lat = Math.max(-85.05112878, Math.min(85.05112878, latDeg));
        double x = R * Math.toRadians(lonDeg);
        double y = R * Math.log(Math.tan(Math.PI / 4 + Math.toRadians(lat) / 2));
        return new double[] {x, y};
    }

    public static void main(String[] args) {
        double[][] pts = { {116.4074, 39.9042}, {0, 0}, {0, 60}, {0, 85} };
        for (double[] p : pts) {
            double[] xy = forward(p[0], p[1]);
            double k = 1.0 / Math.cos(Math.toRadians(p[1]));
            System.out.printf("lat=%.1f  x=%.1f  y=%.1f  scale=%.3f  area=%.3f%n",
                p[1], xy[0], xy[1], k, k * k);
        }
        // 北京附近：scale≈1.305，面积放大≈1.70 倍
    }
}

import org.locationtech.jts.geom.Envelope;
import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.geom.GeometryFactory;
import org.locationtech.jts.geom.util.AffineTransformation;

public class TileDemo {
    static final double R = 6378137.0;

    public static double[] lonLatToMercator(double lon, double lat) {
        double clamped = Math.max(-85.05112878, Math.min(85.05112878, lat));
        double x = R * Math.toRadians(lon);
        double y = R * Math.log(Math.tan(Math.PI / 4 + Math.toRadians(clamped) / 2));
        return new double[] {x, y};
    }

    public static int[] lonLatToTile(double lon, double lat, int z) {
        int n = 1 << z;
        int xt = (int) Math.floor((lon + 180.0) / 360.0 * n);
        double latRad = Math.toRadians(lat);
        int yt = (int) Math.floor(
            (1 - Math.log(Math.tan(latRad) + 1 / Math.cos(latRad)) / Math.PI) / 2 * n);
        return new int[] {xt, yt};
    }

    public static void main(String[] args) {
        double lon = 116.4074;
        double lat = 39.9042;
        int z = 14;
        int[] tile = lonLatToTile(lon, lat, z);
        System.out.println("tile = " + tile[0] + "/" + tile[1] + "/" + z);

        // 瓦片中心点作为局部原点，把几何平移到局部坐标
        double[] origin = lonLatToMercator(lon, lat);
        GeometryFactory gf = new GeometryFactory();
        Geometry local = gf.toGeometry(new Envelope(-1000, 1000, -1000, 1000));
        AffineTransformation toLocal =
            AffineTransformation.translationInstance(-origin[0], -origin[1]);
        Geometry shifted = toLocal.transform(local);
        System.out.println(shifted.getEnvelopeInternal());
        // Envelope 现在以 (0,0) 为中心，数值量级从 1e7 降到 1e3
    }
}

import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.geom.PrecisionModel;
import org.locationtech.jts.io.OrdinateFormat;
import org.locationtech.jts.io.WKTReader;
import org.locationtech.jts.io.WKTWriter;

public class WriteFormatDemo {
    public static void main(String[] args) throws Exception {
        WKTReader r = new WKTReader();
        Geometry g = r.read("POINT (116.40738654 39.90416789)");

        // 1) 默认输出：保留可表示精度
        System.out.println(new WKTWriter().write(g));

        // 2) 按 6 位有效数字精度模型输出
        WKTWriter w1 = new WKTWriter();
        w1.setPrecisionModel(new PrecisionModel(1000000)); // 6 位小数
        System.out.println(w1.write(g));

        // 3) 用 OrdinateFormat 单独格式化坐标值
        OrdinateFormat fmt = OrdinateFormat.create(6);
        System.out.println(fmt.format(116.40738654)); // 116.407387
        System.out.println(fmt.format(39.90416789));  // 39.904168

        // 4) 特殊值有固定文本
        System.out.println(OrdinateFormat.DEFAULT.format(Double.NaN));          // NaN
        System.out.println(OrdinateFormat.DEFAULT.format(Double.POSITIVE_INFINITY)); // Inf
    }
}

import org.locationtech.jts.geom.Envelope;
import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.geom.GeometryCollection;

public final class CrsGuard {
    private CrsGuard() {}

    /** 判断包络是否像经纬度数据（范围 + 非投影量级） */
    public static boolean looksLikeLonLat(Geometry g) {
        if (g.isEmpty()) return true;
        Envelope env = g.getEnvelopeInternal();
        return env.getMinX() >= -180 && env.getMaxX() <= 180
            && env.getMinY() >= -90 && env.getMaxY() <= 90;
    }

    /** 判断是否存在反经线风险：包络宽度过大 */
    public static boolean crossesAntimeridian(Geometry g) {
        return !g.isEmpty() && g.getEnvelopeInternal().getWidth() > 180;
    }

    /** 基础元数据校验：SRID 非 0、坐标有限、类型一致 */
    public static void assertUsable(Geometry g, int expectedSrid) {
        if (g == null) throw new IllegalArgumentException("几何为 null");
        if (g.getSRID() != expectedSrid) {
            throw new IllegalArgumentException(
                "SRID 不匹配: 期望 " + expectedSrid + "，实际 " + g.getSRID());
        }
        if (g instanceof GeometryCollection) {
            for (int i = 0; i < g.getNumGeometries(); i++) {
                assertUsable(g.getGeometryN(i), expectedSrid);
            }
        }
        if (!looksLikeFinite(g)) {
            throw new IllegalArgumentException("坐标含 NaN/Inf");
        }
    }

    private static boolean looksLikeFinite(Geometry g) {
        for (int i = 0; i < g.getNumPoints(); i++) {
            double x = g.getCoordinates()[i].x;
            double y = g.getCoordinates()[i].y;
            if (Double.isNaN(x) || Double.isInfinite(x)) return false;
            if (Double.isNaN(y) || Double.isInfinite(y)) return false;
        }
        return true;
    }
}