等经纬度影像 转 墨卡托影像，等经纬度图片 转 墨卡托图片

前段时间做了气象方面的工作，要在gis地图上叠加卫星影像。

gis地图是墨卡托投影，卫星影像是等经纬度投影，不能直接叠加。

于是整出了这个等经纬度影像 转 墨卡托影像的程序。

上代码：

import java.awt.image.BufferedImage;
import java.io.File;
import java.io.IOException;
import java.util.ArrayList;
import java.util.List;

import javax.imageio.ImageIO;

/**
 * 等经纬度影像 转 墨卡托影像
 * 基本思路：
 *     1、赤道上经度跨越1度，其代表的距离恒定不变。
 *     2、（经纬度影像）由起止经纬度可以计算出径向距离，得到径向有多少像素，从而可以计算出一像素代表多少米；
 *     3、纬度到赤道的距离有公式可以计算出，再由第2步得到的一像素代表多少米，就可以知道纬度需要多少像素填充。
 * @author RainingTime
 *
 */
public class LatLonImage2MercatorImage {
	
	public static void main(String[] args) throws IOException {
		long start = System.currentTimeMillis();
        
		//原始数据起止经纬度{{开始纬度、开始经度}、{结束纬度、结束经度}}，左上角开始，右下角结束
		Double[][] sourceRange = {{53.817,60.2}, {4.83,138.78}};
		
		//裁剪区域起止经纬度{{开始纬度、开始经度}、{结束纬度、结束经度}}，左上角开始，右下角结束。
		//若不需要裁剪，则令cutRange = sourceRange
        Double[][] cutRange = {{53.817,60.2}, {4.83,138.78}};
        
        String imageType = "png";
        
        String source = "C:/Users/RainingTime/Desktop/FY4影像/FY4A_20210421074900000.jpg";
        String outPut = "C:/Users/RainingTime/Desktop/FY4影像/FY4A_20210421074900000-mkt." + imageType;
        BufferedImage sourceImage = ImageIO.read(new File(source));
        BufferedImage mercator = imageToMercator(sourceImage, sourceRange, cutRange, imageType);
		ImageIO.write(mercator, imageType, new File(outPut));
        
        long end = System.currentTimeMillis();
        System.out.println("转换完成，用时："+(end-start)+"ms");
	}

	/**
	 * png、jpg影像转墨卡托投影
	 * @param sourcePath	原始影像
	 * @param outPutPath	输出影像
	 * @param sourceRange	原始影像经纬度范围，要求左上角开始，右下角结束：{{开始纬度、开始经度}、{结束纬度、结束经度}}
	 * @param cutRange		裁剪区域经纬度范围，要求左上角开始，右下角结束：{{开始纬度、开始经度}、{结束纬度、结束经度}}，为null时默认不裁剪
	 * @param imageType		输出影像格式，png、jpg两种
	 * @throws IOException
	 */
	public static BufferedImage imageToMercator(BufferedImage sourceImage, Double[][] sourceRange, Double[][] cutRange, String imageType) throws IOException {
		double sourceStartLat = sourceRange[0][0];
		double sourceStartLon = sourceRange[0][1];
		double sourceEndLat = sourceRange[1][0];
		double sourceEndLon = sourceRange[1][1];
		double cutStartLat = sourceStartLat;
		double cutStartLon = sourceStartLon;
		double cutEndLat = sourceEndLat;
		double cutEndLon = sourceEndLon;
		if(cutRange != null){
			cutStartLat = cutRange[0][0];
			cutStartLon = cutRange[0][1];
			cutEndLat = cutRange[1][0];
			cutEndLon = cutRange[1][1];
		}
		
		if(sourceEndLat>sourceStartLat || sourceStartLon>sourceEndLon){
			System.out.println("数据源经纬度坐标从左上角开始，右下角结束！");
			return null;
		}
		if(cutEndLat>cutStartLat || cutStartLon>cutEndLon){
			System.out.println("裁剪区域经纬度坐标从左上角开始，右下角结束！");
			return null;
		}
		/**
		 * earthRadius = 6378137;//地球赤道半径6378137米
		 * 20037508.3427892 = earthRadius * (math.pi - 0);//赤道周长的一半
		 * 85.05112877980659 = (math.atan(math.exp(aa / earthRadius))-math.pi/4)*2 * 180 / math.pi;//墨卡托最大有效纬度
		 */
		if(cutStartLat > 85.051128){
			System.out.println("墨卡托投影起始纬度最大为：85.051128，裁剪区域将被限定。");
			cutStartLat = 85.051128;
		}
		if(cutEndLat <- 85.051128){
			System.out.println("墨卡托投影终止纬度最小为：-85.051128，裁剪区域将被限定。");
			cutEndLat = -85.051128;
		}
		
        int sourceWidth = sourceImage.getWidth();
		int sourceHeight = sourceImage.getHeight();
		double lonStep = (sourceEndLon-sourceStartLon)/(sourceWidth-1);
		double latStep = (sourceEndLat-sourceStartLat)/(sourceHeight-1);
		System.out.println("sourceWidth: " + sourceWidth + ", sourceHeight: " + sourceHeight + ", latStep: " + latStep + ", lonStep: " + lonStep);

		int cutWidth = (int)((cutEndLon - cutStartLon)/lonStep + 0.5) + 1;
		int cutHeight = (int)((cutEndLat - cutStartLat)/latStep + 0.5) + 1;
		System.out.println("cutWidth: "+cutWidth+", cutHeight: "+cutHeight);
		
		//截止纬度与赤道之间有多少个点，北半球为正值
		int start = (int)((Math.log(Math.tan((90 + cutStartLat) * Math.PI / 360)) / (Math.PI / 180))/ lonStep + 0.5);
		//起始纬度与赤道之间有多少个点，南半球为负值
		int end = (int)((Math.log(Math.tan((90 + cutEndLat) * Math.PI / 360)) / (Math.PI / 180))/ lonStep + 0.5);
		
		int mwidth = cutWidth;//墨卡托投影宽度
		int mheight = start - end + 1;//墨卡托投影高度
		System.out.println("mwidth: "+mwidth+", mheight: "+mheight);

		//建立预设
		List<List<Integer>> list = new ArrayList<List<Integer>>();
		List<Integer> innerList = null;
		for(int y=0; y<mheight; y++){
			innerList = new ArrayList<Integer>();
			for(int x=0; x<mwidth; x++){
				innerList.add(0);//填充0，透明图层
			}
			list.add(innerList);
			innerList = null;
		}
		
		//填充真实数据
		int equator = -1;//赤道的像素y轴坐标
		BufferedImage outImage = null;
		if("png".equals(imageType)){
			outImage = new BufferedImage(mwidth, mheight, BufferedImage.TYPE_INT_ARGB);
		}else{
			outImage = new BufferedImage(mwidth, mheight, BufferedImage.TYPE_INT_RGB);
		}
		List<Integer> yIndexList = new ArrayList<Integer>(); 
		for(int y=0; y<cutHeight; y++){
			double lat = cutStartLat + (y * latStep);//截止纬度 减 当前纬度，得到纬度差
			if(-85.051128 <= lat && lat <= 85.051128){
				int oy = start - (int)((Math.log(Math.tan((90 + lat) * Math.PI / 360)) / (Math.PI / 180)) / lonStep + 0.5);
				
				if(oy<mheight){
					yIndexList.add(oy);
					if(lat == 0.0 || (lat>0.0 && (lat+0.5*latStep)<0.0) || (lat<0.0 && (lat-0.5*latStep)>0.0)){//赤道
						equator = oy;
					}
					for(int x=0; x<cutWidth; x++){
		        		int px = (int)((cutStartLon+x*lonStep-sourceStartLon)/lonStep + 0.5);
		        		int py = (int)((cutStartLat+y*latStep-sourceStartLat)/latStep + 0.5);
						int rgb = sourceImage.getRGB(px, py);
						list.get(oy).set(x, rgb);
					}
				}
			}
		}
		System.out.println("赤道的像素y轴坐标："+equator);
		
		for(int y=0; y<mheight; y++){
			if(y <= equator){
				for(int i=0; i<yIndexList.size(); i++){
					if(y == yIndexList.get(i)){
						for(int x=0; x<cutWidth; x++){
							outImage.setRGB(x, y, list.get(y).get(x));
						}
						break;
					} else if(y > yIndexList.get(i) && (i+1)<yIndexList.size() && y < yIndexList.get(i+1)){
						for(int x=0; x<cutWidth; x++){
							outImage.setRGB(x, y, list.get(yIndexList.get(i+1)).get(x));
						}
						break;
					}
				}
			} else {
				for(int i=0; i<yIndexList.size(); i++){
					if(y == yIndexList.get(i)){
						for(int x=0; x<cutWidth; x++){
							outImage.setRGB(x, y, list.get(y).get(x));
						}
						break;
					} else if(y > yIndexList.get(i) && (i+1)<yIndexList.size() && y < yIndexList.get(i+1)){
						for(int x=0; x<cutWidth; x++){
							outImage.setRGB(x, y, list.get(yIndexList.get(i)).get(x));
						}
						break;
					}
				}
			}
		}
		list = null;
		return outImage;
	}

}

原图：风云四A星真彩色云图，来源：http://www.nmc.cn/publish/satellite/FY4A-true-color.htm

输出：因为不能上传大于5M的图片，这里将图片压缩了下才上传的。

传送门：墨卡托影像 转 等经纬度影像，墨卡托图片 转 等经纬度图片