package cn.service.hspd.utils;import com.alibaba.fastjson2.JSON;import java.awt.geom.Point2D;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.List;
import java.util.Map;
import java.util.stream.Collectors;public class GisUtil {// 两点距离 x,lon 经度 y,lat 纬度private static final double EARTH_RADIUS = 6371000;/*** 计算两点jul** @param p1* @param p2* @return*/public static double jl(Point2D p1, Point2D p2) {double lat1Rad = Math.toRadians(p1.getY());double lon1Rad = Math.toRadians(p1.getX());double lat2Rad = Math.toRadians(p2.getY());double lon2Rad = Math.toRadians(p2.getX());double deltaLat = lat2Rad - lat1Rad;double deltaLon = lon2Rad - lon1Rad;double a = Math.sin(deltaLat / 2) * Math.sin(deltaLat / 2)+ Math.cos(lat1Rad) * Math.cos(lat2Rad)* Math.sin(deltaLon / 2) * Math.sin(deltaLon / 2);double c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1 - a));return EARTH_RADIUS * c;}/*** 判断两条线是否相交返回相交点** @param p1* @param p2* @param p3* @param p4* @return*/public static Point2D.Double xj(Point2D p1, Point2D p2, Point2D p3, Point2D p4) {double x1 = p1.getX(), y1 = p1.getY();double x2 = p2.getX(), y2 = p2.getY();double x3 = p3.getX(), y3 = p3.getY();double x4 = p4.getX(), y4 = p4.getY();double d = (x1 - x2) * (y3 - y4) - (y1 - y2) * (x3 - x4);if (d == 0) return null;double xi = ((x3 - x4) * (x1 * y2 - y1 * x2) - (x1 - x2) * (x3 * y4 - y3 * x4)) / d;double yi = ((y3 - y4) * (x1 * y2 - y1 * x2) - (y1 - y2) * (x3 * y4 - y3 * x4)) / d;if (x1 < Math.min(x1, x2) || xi > Math.max(x1, x2) || xi < Math.min(x3, x4) || xi > Math.max(x3, x4))return null;if (yi < Math.min(y1, y2) || yi > Math.max(y1, y2) || yi < Math.min(y3, y4) || yi > Math.max(y3, y4))return null;return new Point2D.Double(xi, yi);}/*** 中心点** @param list* @return*/public static Point2D.Double zxd(List<Point2D.Double> list) {if (list.size() <= 1) return null;Point2D top = list.get(0);Point2D bot = list.get(1);if (list.size() == 2) {return new Point2D.Double((top.getX() + bot.getX()) / 2.0, (top.getY() + bot.getY()) / 2.0);}double maxJl = -1;for (int i = 0; i < list.size(); i++) {for (int j = i + 1; j < list.size(); j++) {if (Math.abs(i - j) == 1 || (i == 0 && j == list.size() - 1)) continue;double jl = jl(list.get(i), list.get(j));if (jl > maxJl) {maxJl = jl;top = list.get(i);bot = list.get(j);}}}return new Point2D.Double((top.getX() + bot.getX()) / 2.0, (top.getY() + bot.getY()) / 2.0);}/*** 最大直径** @param list* @return*/public static double zdzj(List<Point2D.Double> list) {if (list.size() <= 1) return 0;Point2D top = list.get(0);Point2D bot = list.get(1);if (list.size() == 2) {return jl(top, bot);}double maxJl = -1;for (int i = 0; i < list.size(); i++) {for (int j = i + 1; j < list.size(); j++) {if (Math.abs(i - j) == 1 || (i == 0 && j == list.size() - 1)) continue;double jl = jl(list.get(i), list.get(j));if (jl > maxJl) {maxJl = jl;top = list.get(i);bot = list.get(j);}}}return jl(top, bot);}/*** 最x小直径** @param list* @return*/public static double zxzj(List<Point2D.Double> list) {if (list.size() <= 1) return 0;Point2D top = list.get(0);Point2D bot = list.get(1);if (list.size() == 2) {return jl(top, bot);}double min = jl(list.get(0), list.get(3));for (int i = 0; i < list.size(); i++) {for (int j = i + 1; j < list.size(); j++) {if (Math.abs(i - j) == 1 || (i == 0 && j == list.size() - 1)) continue;double jl = jl(list.get(i), list.get(j));if (jl <= min) {min = jl;top = list.get(i);bot = list.get(j);}}}return jl(top, bot);}public static double zc(List<Point2D.Double> list) {double zc = 0;for (int i = 1; i < list.size(); i++) {zc += jl(list.get(i - 1), list.get(i));}return zc;}/*** 计算体积 圆柱 锥形** @param bj* @param height* @param lx 1 圆柱 2 圆锥* @return*/public static double tj(double bj, double height, int lx) {if (lx == 1) {return Math.PI * Math.pow(bj, 2) * height;}return (1.0 / 3) * Math.PI * Math.pow(bj, 2) * height;}public static double format(Double d, int i) {if (d == null) return 0.0;String format = String.format("%." + i + "f", d);return Double.parseDouble(format);}/*** 像素转经纬坐标** @param zs 图片左上位置* @param yx 图片右下位置* @param width 图片宽* @param height 图片高* @param zb 像素信息* @return*/public static Point2D.Double xsTojw(Point2D zs, Point2D yx, double width, double height, Point2D zb) {double lat = zs.getY() - (zb.getY() * (zs.getY() - yx.getY()) / height);double lon = zs.getX() - (zb.getX() * (yx.getX() - zs.getX()) / width);return new Point2D.Double(lon, lat);}// 长半轴private static final double WGS84_A = 6378137.0;// 短半轴private static final double WGS84_B = 6356752.314245;// 扁率private static final double WGS84_F = 1.0 / 298.257223563;// 第一偏心率平方private static final double WGS84_E2 = 6.69437999014e-3;// 收敛经度private static final double WGS84_EPSILON = 1.0e-12;/*** 三维直角坐标转经纬度** @param x* @param y* @param z* @return*/public static Point2D.Double swTojw(double x, double y, double z) {double longitude = Math.atan2(y, x);double p = Math.sqrt(x * x + y * y);double latitude = Math.atan2(z, p * (1 - WGS84_E2));double height = 0.0;int maxIterations = 10;for (int i = 0; i < maxIterations; i++) {double sinLat = Math.sin(latitude);double n = WGS84_A / Math.sqrt(1 - WGS84_E2 * sinLat * sinLat);double newHeight = p / Math.cos(latitude) - n;double newLatitude = Math.atan2(z, p * (1 - WGS84_E2 * n / (n + newHeight)));if (Math.abs(newLatitude - latitude) < WGS84_EPSILON && Math.abs(newHeight - height) < WGS84_EPSILON) {break;}latitude = newLatitude;height = newHeight;}double lon = Math.toDegrees(longitude);double lat = Math.toDegrees(latitude);return new Point2D.Double(lon, lat);}/*** 经纬高转为三维直角坐标** @param lon 经度* @param lat 纬度* @param height 高度* @return*/public static List<Double> jwTosw(double lon, double lat, double height) {double lonRad = Math.toRadians(lon);double latRad = Math.toRadians(lat);double sinLat = Math.sin(latRad);double coslat = Math.cos(latRad);double sinLon = Math.sin(lonRad);double cosLon = Math.cos(lonRad);double n = WGS84_A / Math.sqrt(1 - WGS84_E2 * sinLat * sinLat);double x = (n + height) * coslat * cosLon;double y = (n + height) * coslat * sinLon;double z = (n * (1 - WGS84_E2) + height) * sinLat;return Arrays.asList(x, y, z);}/*** 坐标经纬转像素** @param zs 图片左上位置* @param yx 图片右下位置* @param width 图片宽* @param height 图片高* @param zb 坐标经纬度* @return*/public static Point2D.Double jwToxs(Point2D zs, Point2D yx, double width, double height, Point2D zb) {int x = (int) (zb.getX() - zs.getX() * width / (yx.getX() - zs.getX()));int y = (int) ((zs.getY() - zb.getY()) * height / (zs.getY() - yx.getY()));return new Point2D.Double(x, y);}/*** 判断当前位置是否在多边形区域内** @param orderLocation 当前点* @param partitionLocation 区域顶点* @return*/public static boolean isInPolygon(Map orderLocation, String partitionLocation) {double p_x = Double.parseDouble((String) orderLocation.get("X"));double p_y = Double.parseDouble((String) orderLocation.get("Y"));Point2D.Double point = new Point2D.Double(p_x, p_y);List<Point2D.Double> pointList = new ArrayList<Point2D.Double>();List objects = JSON.parseArray(partitionLocation);for (Object str : objects) {List<Double> doubles = JSON.parseArray(JSON.toJSONString(str), Double.class);double polygonPoint_x = doubles.get(0);double polygonPoint_y = doubles.get(1);Point2D.Double polygonPoint = new Point2D.Double(polygonPoint_x, polygonPoint_y);pointList.add(polygonPoint);}return IsPtInPoly(point, pointList);}/*** 返回一个点是否在一个多边形区域内， 如果点位于多边形的顶点或边上，不算做点在多边形内，返回false** @param point* @param polygon* @return*/public static boolean checkWithJdkGeneralPath(Point2D.Double point, List<Point2D.Double> polygon) {java.awt.geom.GeneralPath p = new java.awt.geom.GeneralPath();Point2D.Double first = polygon.get(0);p.moveTo(first.x, first.y);polygon.remove(0);for (Point2D.Double d : polygon) {p.lineTo(d.x, d.y);}p.lineTo(first.x, first.y);p.closePath();return p.contains(point);}/*** 判断点是否在多边形内，如果点位于多边形的顶点或边上，也算做点在多边形内，直接返回true** @param point 检测点* @param pts 多边形的顶点* @return 点在多边形内返回true, 否则返回false*/public static boolean IsPtInPoly(Point2D.Double point, List<Point2D.Double> pts) {int N = pts.size();boolean boundOrVertex = true; //如果点位于多边形的顶点或边上，也算做点在多边形内，直接返回trueint intersectCount = 0;//cross points count of xdouble precision = 2e-10; //浮点类型计算时候与0比较时候的容差Point2D.Double p1, p2;//neighbour bound verticesPoint2D.Double p = point; //当前点p1 = pts.get(0);//left vertexfor (int i = 1; i <= N; ++i) {//check all raysif (p.equals(p1)) {return boundOrVertex;//p is an vertex}p2 = pts.get(i % N);//right vertexif (p.x < Math.min(p1.x, p2.x) || p.x > Math.max(p1.x, p2.x)) {//ray is outside of our interestsp1 = p2;continue;//next ray left point}if (p.x > Math.min(p1.x, p2.x) && p.x < Math.max(p1.x, p2.x)) {//ray is crossing over by the algorithm (common part of)if (p.y <= Math.max(p1.y, p2.y)) {//x is before of rayif (p1.x == p2.x && p.y >= Math.min(p1.y, p2.y)) {//overlies on a horizontal rayreturn boundOrVertex;}if (p1.y == p2.y) {//ray is verticalif (p1.y == p.y) {//overlies on a vertical rayreturn boundOrVertex;} else {//before ray++intersectCount;}} else {//cross point on the left sidedouble xinters = (p.x - p1.x) * (p2.y - p1.y) / (p2.x - p1.x) + p1.y;//cross point of yif (Math.abs(p.y - xinters) < precision) {//overlies on a rayreturn boundOrVertex;}if (p.y < xinters) {//before ray++intersectCount;}}}} else {//special case when ray is crossing through the vertexif (p.x == p2.x && p.y <= p2.y) {//p crossing over p2Point2D.Double p3 = pts.get((i + 1) % N); //next vertexif (p.x >= Math.min(p1.x, p3.x) && p.x <= Math.max(p1.x, p3.x)) {//p.x lies between p1.x & p3.x++intersectCount;} else {intersectCount += 2;}}}p1 = p2;//next ray left point}if (intersectCount % 2 == 0) {//偶数在多边形外return false;} else { //奇数在多边形内return true;}}
}
package cn.service.hspd.utils;import java.awt.geom.AffineTransform;
import java.awt.geom.Point2D;public class GeoPixelConverter {private final AffineTransform transform;private final double pixelError;public GeoPixelConverter(double topLeftLon, double topLeftLat,double bottomRightLon, double bottomRightLat,int widthPixels, int heightPixels) {// 计算仿射变换参数（保留10位小数）double scaleX = (bottomRightLon - topLeftLon) / widthPixels;double scaleY = (bottomRightLat - topLeftLat) / heightPixels;this.transform = new AffineTransform(scaleX, 0, 0, scaleY, topLeftLon, topLeftLat);// 理论最大误差计算（基于双精度浮点运算）this.pixelError = Math.max(Math.abs(scaleX) * 1e-10,Math.abs(scaleY) * 1e-10);}public Point2D.Double pixelToGeo(int x, int y) {Point2D.Double src = new Point2D.Double(x, y);Point2D.Double dst = new Point2D.Double();transform.transform(src, dst);// 结果四舍五入到小数点后8位dst.x = Math.round(dst.x * 1e8) / 1e8;dst.y = Math.round(dst.y * 1e8) / 1e8;return dst;}public Point2D geoToPixel(double longitude, double latitude) {try {Point2D.Double src = new Point2D.Double(longitude, latitude);Point2D.Double dst = new Point2D.Double();transform.inverseTransform(src, dst);// 像素坐标取整dst.x = Math.round(dst.x);dst.y = Math.round(dst.y);return dst;} catch (Exception e) {throw new IllegalArgumentException("坐标超出转换范围");}}public double getMaxError() {return pixelError;}// 边缘像素补偿算法示例public static Point2D.Double adjustEdgePixel(int x, int y, int width, int height) {// 对边缘像素进行0.5像素偏移补偿double adjX = (x == 0) ? 0.5 : (x == width-1) ? width-1.5 : x;double adjY = (y == 0) ? 0.5 : (y == height-1) ? height-1.5 : y;return new Point2D.Double(adjX, adjY);}public static void main(String[] args) {// 示例：故宫影像图转换（误差<0.00001度）GeoPixelConverter converter = new GeoPixelConverter(116.391111, 39.916667, // 左上角经纬度116.403889, 39.908333, // 右下角经纬度4096, 4096); // 图像宽高Point2D.Double geo = converter.pixelToGeo(4096, 0);System.out.printf("中心点坐标: %.8f, %.8f%n", geo.x, geo.y);System.out.printf("理论最大误差: %.10f 度", converter.getMaxError());}
}
