总体思路：
1.统一输入坐标数据，如果是 4326 转为 3857，不管是单个 还是一组，统一转为 一组
2.计算瓦片坐标XY，基于3857 坐标（Web墨卡托或Web Mercator），获取 22 级的瓦片 X、Y，（为什么用22级，因为我们只生成了 22级的瓦片数据）
3.获取坐标相对瓦片的位置，计算坐标在单张瓦片图上的的相对位置 RelativeX RelativeY，左上为0,0
4.获取灰度值，根据瓦片XY 以及 相对位置 RelativeX RelativeY，获取坐标对应瓦片图的像素位置上的像素值
5.获取单张瓦片基准值
6.获取高程，利用像素值 和 瓦片基准值 进行高程的计算
背景说明
我们只提供 22 级的瓦片图
package com.sxkc.modules.util;
import com.alibaba.fastjson.JSONObject;
import com.google.common.collect.Lists;
import lombok.extern.slf4j.Slf4j;
import org.apache.commons.lang.StringUtils;
import org.apache.commons.lang3.ObjectUtils;
import org.geotools.referencing.CRS;
import org.locationtech.jts.geom.Coordinate;
import org.locationtech.jts.geom.GeometryFactory;
import org.opengis.referencing.crs.CRSAuthorityFactory;
import org.opengis.referencing.crs.CoordinateReferenceSystem;
import org.opengis.referencing.operation.MathTransform;
import javax.imageio.ImageIO;
import java.awt.image.*;
import java.io.File;
import java.io.FileReader;
import java.io.Reader;
import java.math.BigDecimal;
import java.math.RoundingMode;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.List;
import java.util.concurrent.Callable;
import java.util.concurrent.ExecutorService;
import java.util.concurrent.Executors;
import java.util.concurrent.Future;
import java.util.stream.Collectors;
import com.google.gson.Gson;
import org.springframework.beans.factory.annotation.Value;
import org.springframework.stereotype.Component;
@Slf4j
@Component
public class HeightUtil {
@Value("${file.resourcePath}")
private String resourcePath;
// private static final String JSON_FILEs = "%s/elv/1/22/%s/%s.json";
// private static final String TILE_FILEs = "%s/elv/1/22/%s/%s.png";
private static final String JSON_FILE1 = "%s/elv/";
private static final String TILE_FILE1 = "%s/elv/";
private static final String JSON_FILE2 = "/22/%s/%s.json";
private static final String TILE_FILE2 = "/22/%s/%s.png";
private static GeometryFactory geometryFactory = new GeometryFactory();// jts GeometryFactory
private static CoordinateReferenceSystem crs4326; // WGS84 经纬度坐标系
private static CoordinateReferenceSystem crs3857; // 定义 Web Mercator 投影坐标系
private String taskId = "";
private int zoomLevel = 22;
private String inputEpsg = "4326"; // 默认为4326 可传3857
private String JSON_FILE="";
private String TILE_FILE="";
public HeightUtil() {
}
public void init(String taskId) {
this.taskId = taskId;
}
public void init(String taskId, String inputEpsg) {
this.taskId = taskId;
this.inputEpsg = inputEpsg;
}
public void init(String taskId, String inputEpsg,String zlevel) {
this.taskId = taskId;
this.inputEpsg = inputEpsg;
if(StringUtils.isEmpty(zlevel)){
zlevel="0";
}
this.JSON_FILE=JSON_FILE1+zlevel+JSON_FILE2;
this.TILE_FILE=TILE_FILE1+zlevel+TILE_FILE2;
}
// 坐标转换
private double[] lonlat2xy(double[] lonlat){
double[] xy = new double[2];
try {
// 创建坐标系对象
if(ObjectUtils.isEmpty(crs4326) || ObjectUtils.isEmpty(crs3857)){
CRSAuthorityFactory factory = CRS.getAuthorityFactory(true);
crs4326 = factory.createCoordinateReferenceSystem("EPSG:4326");
// 定义 Web Mercator 投影坐标系
crs3857 = CRS.decode("EPSG:3857");
}
// 创建转换器，从 WGS 84 到 Web Mercator
MathTransform transform = CRS.findMathTransform(crs4326, crs3857);
// 将经纬度坐标转换为投影坐标
double[] projectedCoords = new double[2];
transform.transform(new double[]{lonlat[0], lonlat[1]}, 0, projectedCoords, 0, 1);
xy = projectedCoords;
// System.out.println("Web Mercator XY: " + xy[0] + ", " + xy[1]);
}catch (Exception e){
e.printStackTrace();
}
return xy;
}
private Integer[] getTileXY(double[] xy){
Integer[] tileXY = new Integer[2];
try {
// 计算瓦片坐标
int tileX = (int) Math.floor((xy[0] + 20037508.34) / (2 * 20037508.34) * (1 << zoomLevel));
int tileY = (int) Math.floor((20037508.34 - xy[1]) / (2 * 20037508.34) * (1 << zoomLevel));
// System.out.println("Tile Coordinates: " + tileX + ", " + tileY);
tileXY[0] = tileX;
tileXY[1] = tileY;
}catch (Exception e){
e.printStackTrace();
}
return tileXY;
}
// 获取经纬度相对于单个瓦片的相对坐标
private Integer[] getRelativeXY(double[] xy, Integer[] tileXY){
Integer[] relativeXY = new Integer[2];
BigDecimal mercatorX = BigDecimal.valueOf(xy[0]);
BigDecimal mercatorY = BigDecimal.valueOf(xy[1]);
// 计算瓦片的空间范围 参考 mercantile.xy_bounds 的公式 , 使用 BigDecimal 计算更精确
// double CE = 2 * Math.PI * 6378137.0;
// double tileSize = CE / Math.pow(2, zoomLevel);
// double CE_2 = CE/2;
BigDecimal CE = BigDecimal.valueOf(2 * Math.PI * 6378137.0);
BigDecimal tileSize = CE.divide(BigDecimal.valueOf(Math.pow(2, zoomLevel)), 13 , BigDecimal.ROUND_UP);
BigDecimal CE_2 = CE.divide(BigDecimal.valueOf(2), 13 , BigDecimal.ROUND_UP);
BigDecimal minX = tileSize.multiply(BigDecimal.valueOf(tileXY[0])).subtract(CE_2);
BigDecimal maxX = minX.add(tileSize);
BigDecimal maxY = CE_2.subtract(tileSize.multiply(BigDecimal.valueOf(tileXY[1])));
BigDecimal minY = maxY.subtract(tileSize);
// System.out.println("Tile Spatial Range: " + minX + ", " + minY + ", " + maxX + ", " + maxY);
// 利用 （目标点X - 最小X）/ （最大X - 最小X）计算出x比例，乘以 256 就得到了基于图片的坐标
BigDecimal ratioX = (mercatorX.subtract(minX)).divide(maxX.subtract(minX), 10, RoundingMode.HALF_UP );
BigDecimal ratioY = (BigDecimal.valueOf(1).subtract((mercatorY.subtract(minY)).divide(maxY.subtract(minY), 10, RoundingMode.HALF_UP )));
BigDecimal relativeX = ratioX.multiply(BigDecimal.valueOf(256));
BigDecimal relativeY = ratioY.multiply(BigDecimal.valueOf(256));
// System.out.println("Tile relative ratio: " + ratioX + ", " + ratioY );
// System.out.println("Tile relative XY: " + relativeX + ", " + relativeY );
relativeXY[0] = (int) relativeX.floatValue();
relativeXY[1] = (int) relativeY.floatValue();
return relativeXY;
}
// 获取灰度值
private Integer getGray(Integer[] tileXY, Integer relativeXY[]){
String tileImagePath = String.format(resourcePath + TILE_FILE, this.taskId, tileXY[0], tileXY[1]);
Integer gray = 0;
try {
// 读取图像文件
File imageFile = new File(tileImagePath);
BufferedImage image = ImageIO.read(imageFile);
// 获取单通道图像的指定像素位置的灰度值
Raster raster = image.getRaster();
gray = raster.getSample(relativeXY[0], relativeXY[1], 0);
// 输出像素值
// System.out.println("灰度值：" + gray);
return gray;
} catch (Exception e) {
// e.printStackTrace();
}
return gray;
}
// 获取elv
private Float getElv(Integer[] tileXY){
Float elv = 0F;
String tileJsonPath = String.format(resourcePath + JSON_FILE, this.taskId, tileXY[0], tileXY[1]);
try {
Gson gson = new Gson();
Reader reader = new FileReader(tileJsonPath);
JSONObject object = gson.fromJson(reader, JSONObject.class);
if (object.containsKey("elv")){
elv = Float.valueOf(object.get("elv").toString());
}
// 输出elv
// System.out.println("elv：" + elv);
} catch (Exception e) {
// e.printStackTrace();
}
return elv;
}
public Float getHeight(double longitude, double latitude){
Float height = 0f;
try {
double[] lonlat = null;
double[] xy = null;
if(this.inputEpsg.equals("4326")) {
lonlat = new double[]{longitude, latitude};
// 墨卡托坐标
xy = lonlat2xy(lonlat);
}else{
xy = new double[]{longitude, latitude};
}
// 通过 墨卡托坐标 过去瓦片 xy
Integer[] tileXY = getTileXY(xy);
// 获取经纬度 相对于瓦片的相对坐标 , 从而获取灰度值
Integer[] relativeXY = getRelativeXY(xy, tileXY);
Integer gray = getGray(tileXY, relativeXY);
// 获取elv
Float elv = getElv(tileXY);
height = (gray/(255/60)-30)/100 + elv;
}catch (Exception e) {
// e.printStackTrace();
}
return height;
}
public Coordinate[] getHeight(Coordinate[] coordinates){
List<double[]> coordinateList = Arrays.stream(coordinates)
.map(coordinate -> new double[]{coordinate.getX(), coordinate.getY()})
.collect(Collectors.toList());
List<BigDecimal[]> coordinateListZ = this.getHeight(coordinateList);
Coordinate[] coordinatesReturn = coordinateListZ.stream()
.map(bigDecimalArray -> new Coordinate(bigDecimalArray[0].doubleValue(), bigDecimalArray[1].doubleValue(), bigDecimalArray[2].doubleValue()))
.toArray(Coordinate[]::new);
return coordinatesReturn;
}
public List<BigDecimal[]> getHeight(List<double[]> inputCoords){
List<BigDecimal[]> coordinates = new ArrayList<>();
// for (double[] item : inputCoords) {
// float z = getHeight(item[0], item[1]);
// coordinates.add(new double[]{item[0], item[1], z});
// }
// if(coordinates.size() > 0){
// return coordinates;
// }
// 每组数量，最多分10组
Integer groupNum = (int)(inputCoords.size()/10) + 1;
List<List<double[]>> inputCoordsParttion = Lists.partition(inputCoords, groupNum);
try {
// 创建线程池
ExecutorService executorService = Executors.newFixedThreadPool(inputCoordsParttion.size());
// 创建任务列表
List<Callable<List<BigDecimal[]>>> tasks = new ArrayList<>();
Integer i = 0;
for (List<double[]> item : inputCoordsParttion) {
tasks.add(new DataProcessingTask(this, item));
i++;
}
// 提交任务并获取Future对象列表
List<Future<List<BigDecimal[]>>> futures = executorService.invokeAll(tasks);
// 处理任务结果 这里是阻塞式的 有先后顺序的返回
for (Future<List<BigDecimal[]>> future : futures) {
// 获取任务结果，注意此处会阻塞等待任务完成
List<BigDecimal[]> result = future.get();
coordinates.addAll(result);
}
// 关闭线程池
executorService.shutdown();
}catch (Exception e){
e.printStackTrace();
}
return coordinates;
}
}
class DataProcessingTask implements Callable<List<BigDecimal[]>> {
private List<double[]> coords;
private HeightUtil heightUtil;
public DataProcessingTask(HeightUtil heightUtil, List<double[]> coords) {
this.coords = coords;
this.heightUtil = heightUtil;
}
@Override
public List<BigDecimal[]> call() throws Exception {
List<BigDecimal[]> coordinates = new ArrayList<>();
for (double[] item : coords) {
float z = heightUtil.getHeight(item[0], item[1]);
coordinates.add(new BigDecimal[]{BigDecimal.valueOf(item[0]), BigDecimal.valueOf(item[1]), BigDecimal.valueOf(z)});
}
return coordinates;
}
}