问题遇到的现象和发生背景
目的：实现地图瓦片加偏
数据源：EPSG:4326 .TIF 输出：EPSG:3857 .PNG
目前是有两种策略：
1.是直接对数据源进行加偏，然后进行切片（这种方式以验证有效）
@Test
void convertData() {
Dataset srcDs = gdal.Open("D:\\Users\\huaxiaoyang.li\\Desktop\\Image\\19_1.TIF");
Driver TIFDiver = gdal.GetDriverByName("GTIFF");
String s = srcDs.GetProjection();
double[] doubles = srcDs.GetGeoTransform();
// WGS84ToGCJ_02加偏函数已验证准确
double[] doubles1 = CoordinateTransformUtil.WGS84ToGCJ_02(doubles[0], doubles[3]);
Dataset outds = TIFDiver.CreateCopy("D:\\Users\\huaxiaoyang.li\\Desktop\\out.tif", srcDs);
outds.SetProjection(s);
outds.SetGeoTransform(new double[]{doubles1[0], doubles[1], doubles[2], doubles1[1], doubles[4], doubles[5]});
srcDs.delete();
outds.delete();
}
public class CoordinateTransformUtil {
protected final static double SemiMajorAxis = 6378137.0;
protected final static double EarthCircumference = 2 * Math.PI * SemiMajorAxis;
protected final static double mercatorMax = EarthCircumference / 2;
private static final double PI = 3.141592653589793;
private static final double a = 6378245.0;
private static final double f = 1 / 298.3;
private static final double b = a * (1 - f);
private static final double ee = 1 - (b * b) / (a * a);
/**
* @Author: LiHua
* @Description: check if the point in china
* @DateTime: 2022/8/24 11:33
* @Params:
* @Return
*/
private static boolean outOfChina(double lon, double lat) {
if (72.004 <= lon && lon <= 137.8347 && 0.8293 <= lat && lat <= 55.8271) {
return false;
} else {
return true;
}
}
public static double transformLat(double x, double y) {
double ret = -100.0 + 2.0 * x + 3.0 * y + 0.2 * y * y + 0.1 * x * y + 0.2 * Math.sqrt(Math.abs(x));
ret = ret + (20.0 * Math.sin(6.0 * x * PI) + 20.0 * Math.sin(2.0 * x * PI)) * 2.0 / 3.0;
ret = ret + (20.0 * Math.sin(y * PI) + 40.0 * Math.sin(y / 3.0 * PI)) * 2.0 / 3.0;
ret = ret + (160.0 * Math.sin(y / 12.0 * PI) + 320.0 * Math.sin(y * PI / 30.0)) * 2.0 / 3.0;
return ret;
}
public static double transformLon(double x, double y) {
double ret = 300.0 + x + 2.0 * y + 0.1 * x * x + 0.1 * x * y + 0.1 * Math.sqrt(Math.abs(x));
ret = ret + (20.0 * Math.sin(6.0 * x * PI) + 20.0 * Math.sin(2.0 * x * PI)) * 2.0 / 3.0;
ret = ret + (20.0 * Math.sin(x * PI) + 40.0 * Math.sin(x / 3.0 * PI)) * 2.0 / 3.0;
ret = ret + (150.0 * Math.sin(x / 12.0 * PI) + 300.0 * Math.sin(x * PI / 30.0)) * 2.0 / 3.0;
return ret;
}
/**
* @Author: LiHua
* @Description: WGS84转GCJ_02
* @DateTime: 2022/8/24 12:43
* @Params:
* @Return
*/
public static double[] WGS84ToGCJ_02(double wgsLon, double wgsLat) {
if (outOfChina(wgsLon, wgsLat)) {
return new double[]{wgsLon, wgsLat};
}
double dLat = transformLat(wgsLon - 105.0, wgsLat - 35.0);
double dLon = transformLon(wgsLon - 105.0, wgsLat - 35.0);
double radLat = wgsLat / 180.0 * PI;
double magic = Math.sin(radLat);
magic = 1 - ee * magic * magic;
double sqrtMagic = Math.sqrt(magic);
dLat = (dLat * 180.0) / ((a * (1 - ee)) / (magic * sqrtMagic) * PI);
dLon = (dLon * 180.0) / (a / sqrtMagic * Math.cos(radLat) * PI);
double gcjLat = wgsLat + dLat;
double gcjLon = wgsLon + dLon;
return new double[]{gcjLon, gcjLat};
}
/**
* @Author: LiHua
* @Description: GCJ_02转WGS84坐标
* @DateTime: 2022/8/24 16:52
* @Params:
* @Return
*/
public static double[] GCJ_02ToWGS84(double gcjLon, double gcjLat) {
/*
* todo 目前反算采用Fixed Point Iteration；当对精度或是计算效率有特殊需求时可采用其他计算方法，如：API法、二分法等
* */
double[] g0 = new double[]{gcjLon, gcjLat};
double[] g1 = WGS84ToGCJ_02(g0[0], g0[1]);
double[] w0 = new double[]{gcjLon, gcjLat};
double[] w1 = new double[2];
double[] delta = new double[]{Double.MAX_VALUE, Double.MAX_VALUE};
for (int i = 0; i < g0.length; i++) {
w1[i] = w0[i] - (g1[i] - g0[i]);
}
while (Math.abs(delta[0]) >= 1e-9 || Math.abs(delta[1]) >= 1e-9) {
w0 = new double[]{w1[0], w1[1]};
w1 = new double[2];
delta = new double[2];
g1 = WGS84ToGCJ_02(w0[0], w0[1]);
for (int i = 0; i < g0.length; i++) {
w1[i] = w0[i] - (g1[i] - g0[i]);
delta[i] = w1[i] - w0[i];
}
}
return w1;
}
/**
* @Author: LiHua
* @Description: GCJ_02转百度坐标
* @DateTime: 2022/8/24 17:06
* @Params:
* @Return
*/
public static double[] GCJ_02ToBD(double gcjLon, double gcjLat) {
double z = Math.sqrt(gcjLon * gcjLon + gcjLat * gcjLat) + 0.00002 * Math.sin(gcjLat * PI * 3000.0 / 180.0);
double theta = Math.atan2(gcjLat, gcjLon) + 0.000003 * Math.cos(gcjLon * PI * 3000.0 / 180.0);
double bdLon = z * Math.cos(theta) + 0.0065;
double bdLat = z * Math.sin(theta) + 0.006;
return new double[]{bdLon, bdLat};
}
/**
* @Author: LiHua
* @Description: 百度转GCJ_02坐标
* @DateTime: 2022/8/24 17:07
* @Params:
* @Return
*/
public static double[] BDToGCJ_02(double bdLon, double bdLat) {
double x = bdLon - 0.0065;
double y = bdLat - 0.006;
double z = Math.sqrt(x * x + y * y) - 0.00002 * Math.sin(y * PI * 3000.0 / 180.0);
double theta = Math.atan2(y, x) - 0.000003 * Math.cos(x * PI * 3000.0 / 180.0);
double gcjLon = z * Math.cos(theta);
double gcjLat = z * Math.sin(theta);
return new double[]{gcjLon, gcjLat};
}
/**
* @Author: LiHua
* @Description: WGS84坐标转百度坐标
* @DateTime: 2022/8/24 17:21
* @Params:
* @Return
*/
public static double[] WGS84ToBD(double wgsLon, double wgsLat) {
double[] gcj = WGS84ToGCJ_02(wgsLon, wgsLat);
return GCJ_02ToBD(gcj[0], gcj[1]);
}
/**
* @Author: LiHua
* @Description: 百度坐标转WGS84坐标
* @DateTime: 2022/8/24 17:21
* @Params:
* @Return
*/
public static double[] BDToWGS84(double bdLon, double bdLat) {
double[] gcj = BDToGCJ_02(bdLon, bdLat);
return BDToWGS84(gcj[0], gcj[1]);
}
/**
* @Author: LiHua
* @Description: WebMercator to GCJ_02
* @DateTime: 2022/8/26 9:45
* @Params: WebMercator坐标
* @Return WebMercator火星偏转坐标
*/
public static double[] WMToGCJ_02(double mercatorx, double mercatory) {
double[] doubles = transformTo4326(mercatorx, mercatory);
double[] doubles1 = WGS84ToGCJ_02(doubles[0], doubles[1]);
return transformTo3857(doubles1[0], doubles1[1]);
// return WGS84ToGCJ_02(doubles[0], doubles[1]);
}
/**
* @Author: LiHua
* @Description: GCJ_02 to WebMercator
* @DateTime: 2022/8/26 9:49
* @Params: WebMercator火星偏转坐标
* @Return
*/
public static double[] GCJ_02ToWM(double mercatorx, double mercatory) {
double[] doubles = transformTo4326(mercatorx, mercatory);
double[] doubles1 = GCJ_02ToWGS84(doubles[0], doubles[1]);
return transformTo3857(doubles1[0], doubles1[1]);
}
public final static double[] transformTo4326(double mercatorx, double mercatory) {
double lon = mercatorx / mercatorMax * 180;
double lat = mercatory / mercatorMax * 180;
lat = (180 / PI) * (2 * Math.atan(Math.exp((lat * PI) / 180)) - PI / 2);
return new double[]{lon, lat};
}
public final static double[] transformTo3857(double lon, double lat) {
double mercatorx = lon * mercatorMax / 180;
double temp = lat * PI / 180;
double mercatory = SemiMajorAxis / 2 * Math.log((1.0 + Math.sin(temp)) / (1.0 - Math.sin(temp)));
return new double[]{mercatorx, mercatory};
}
// public static void main(String[] args) {
// double[] doubles1 = WGS84ToGCJ_02(120, 25);
// double[] doubles3 = WGS84ToGCJ_02(125, 25);
// double[] doubles4 = WGS84ToGCJ_02(125, 20);
// double[] doubles2 = WGS84ToGCJ_02(120, 20);
//// double[] doubles1 = WGS84ToGCJ_02(120, 25);
//
// int a = 100;
// }
}
2.是先计算地图瓦片坐标及对应的地理坐标（EPSG3857），然后对瓦片坐标进行偏转，最后对偏转的范围进行切片。
def convertToGCJ_02(inrdd: RDD[DataObjInterface], infile: String): RDD[DataObjInterface] = {
inrdd.mapPartitions((part => {
part.map(p => {
val gridtype = p.getAsTxtByKey("gridtype")
val gridname = p.getAsTxtByKey("gridname")
val minx = p.getAsTxtByKey("minx").toDouble
val miny = p.getAsTxtByKey("miny").toDouble
val maxx = p.getAsTxtByKey("maxx").toDouble
val maxy = p.getAsTxtByKey("maxy").toDouble
//瓦片坐标
val zoomstr = p.getAsTxtByKey("zoom").toInt
val x = p.getAsTxtByKey("x").toLong
val y = p.getAsTxtByKey("y").toLong
val double1: Array[Double] = CoordinateTransformUtil.GCJ_02ToWM(minx, miny)
val double2: Array[Double] = CoordinateTransformUtil.GCJ_02ToWM(maxx, maxy)
val double3: Array[Double] = CoordinateTransformUtil.GCJ_02ToWM(minx, maxy)
val double4: Array[Double] = CoordinateTransformUtil.GCJ_02ToWM(maxx, miny)
val xList = List(double1(0), double2(0), double3(0), double4(0))
val yList = List(double1(1), double2(1), double3(1), double4(1))
val cminx = xList.reduce((d1, d2) => if (d1 > d2) d2 else d1)
val cmaxx = xList.reduce((d1, d2) => if (d1 < d2) d2 else d1)
val cminy = yList.reduce((d1, d2) => if (d1 > d2) d2 else d1)
val cmaxy = yList.reduce((d1, d2) => if (d1 < d2) d2 else d1)
//val grid = GridFactory.createGrid(gridname, 256)
//val extent = grid.getTileExtent(zoomstr, new CRSBounds(cminx, cmaxx, cminy, cmaxy))
// val tiles = grid.getTileExtent(zoomstr, new CRSBounds(minx, miny,maxx, maxy))
val tileobj = new RasterTileDataObj(gridtype, gridname,
zoomstr, x, y,
cminx, cmaxx, cminy, cmaxy)
tileobj.setDatafile(infile)
tileobj
})
}))
}
但是切片后的服务展示效果不能满足预期
与高德底图有明显的偏差
有没有人帮忙分析一下出现这种情况的原因可能是什么？
我想要达到的结果
