结合Geotools实现百度09，国测局02和经纬度的相互转换

概述

本文讲述在Java中，结合结合Geotools实现百度09，国测局02和经纬度shp数据的相互转换。

结果

说明：

1、红色的线条是百度09的；
2、蓝色的线条是国测局02的；
3、填充的是原始wgs84的。
4、从图中可以看出，gcj02和wgs84的区别不是很大在一些不是很精确地情况下可以认为是一样的，bd09的区别稍微大一点；

# 实现思路
由于坐标转换是单个点的，所以在处理一个shp的坐标转换的时候，也是一个个点去做转换的。

实现代码

1.ProjTransform.java

package com.lzugis.geotools.utils;

/**
 * @author lzugis
 * 提供了百度坐标（BD09）、国测局坐标（火星坐标，GCJ02）、和WGS84坐标系之间的转换
 * 命名规则：
 *    1、bd代表百度的坐标，gcj代表国测局火星坐标，wgs代表wgs84坐标
 */
public class ProjTransform {
    /**
     * 定义一些常量
     */
    private final double x_PI = 3.14159265358979324 * 3000.0 / 180.0;
    private final double PI = 3.1415926535897932384626;
    private final double a = 6378245.0;
    private final double ee = 0.00669342162296594323;

    /**
     * 百度坐标系 (BD-09) 与 火星坐标系 (GCJ-02)的转换
     * 即 百度 转 谷歌、高德
     * @param bd_lon
     * @param bd_lat
     * @returns {*[]}
     */
    public double[] bd09togcj02(double bd_lon, double bd_lat){
        double x = bd_lon - 0.0065;
        double y = bd_lat - 0.006;
        double z = Math.sqrt(x * x + y * y) - 0.00002 * Math.sin(y * x_PI);
        double theta = Math.atan2(y, x) - 0.000003 * Math.cos(x * x_PI);
        double gg_lon = z * Math.cos(theta);
        double gg_lat = z * Math.sin(theta);
        return new double[]{gg_lon, gg_lat};
    }
    /**
     * 火星坐标系 (GCJ-02) 与百度坐标系 (BD-09) 的转换
     * 即谷歌、高德 转 百度
     * @param gcj_lon
     * @param gcj_lat
     * @returns {*[]}
     */
    public double[] gcj02tobd09(double gcj_lon, double gcj_lat){
        double z = Math.sqrt(gcj_lon * gcj_lon + gcj_lat * gcj_lat) + 0.00002 * Math.sin(gcj_lat * x_PI);
        double theta = Math.atan2(gcj_lat, gcj_lon) + 0.000003 * Math.cos(gcj_lon * x_PI);
        double bd_lon = z * Math.cos(theta) + 0.0065;
        double bd_lat = z * Math.sin(theta) + 0.006;
        return new double[]{bd_lon, bd_lat};
    }
    /**
     * WGS84转GCj02
     * @param wgs_lon
     * @param wgs_lat
     * @returns {*[]}
     */
    public double[] wgs84togcj02(double wgs_lon, double wgs_lat){
        if (out_of_china(wgs_lon, wgs_lat)) {
            return new double[]{wgs_lon, wgs_lat};
        } else {
            double dlat = transformlat(wgs_lon - 105.0, wgs_lat - 35.0);
            double dlon = transformlon(wgs_lon - 105.0, wgs_lat - 35.0);
            double radlat = wgs_lat / 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 mglat = wgs_lat + dlat;
            double mglon = wgs_lon + dlon;
            return new double[]{mglon, mglat};
        }
    }
    /**
     * GCJ02 转换为 WGS84
     * @param gcj_lon
     * @param gcj_lat
     * @returns {*[]}
     */
    public double[] gcj02towgs84(double gcj_lon, double gcj_lat){
        if (out_of_china(gcj_lon, gcj_lat)) {
            return new double[]{gcj_lon, gcj_lat};
        } else {
            double dlat = transformlat(gcj_lon - 105.0, gcj_lat - 35.0);
            double dlng = transformlon(gcj_lon - 105.0, gcj_lat - 35.0);
            double radlat = gcj_lat / 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);
            dlng = (dlng * 180.0) / (a / sqrtmagic * Math.cos(radlat) * PI);
            double mglat = gcj_lat + dlat;
            double mglng = gcj_lon + dlng;
            return new double[]{gcj_lon * 2 - mglng, gcj_lat * 2 - mglat};
        }
    }
    /**
     * 百度转换为wgs84
     * @param bd_lon
     * @param bd_lat
     * @return
     */
    public double[] bd09towgs84(double bd_lon, double bd_lat){
        //1、bd09->gcj
        double[] bd09_gcj02 = bd09togcj02(bd_lon, bd_lat);
        //2、gcj->wgs84
        return gcj02towgs84(bd09_gcj02[0], bd09_gcj02[1]);
    }
    /**
     * wgs84z转换为百度坐标
     * @param wgs_lon
     * @param wgs_lat
     * @return
     */
    public double[] wgs84tobd09(double wgs_lon, double wgs_lat){
        //1、wgs84->gcj
        double[] wgs84_gcj02 = wgs84togcj02(wgs_lon, wgs_lat);
        //2、gcj->bd09
        return gcj02tobd09(wgs84_gcj02[0], wgs84_gcj02[1]);
    }

    /**
     * 判断是否中国境内
     * @param lon
     * @param lat
     * @return
     */
    private boolean out_of_china(double lon, double lat){
        // 纬度3.86~53.55,经度73.66~135.05
        return !(lon > 73.66 && lon < 135.05 && lat > 3.86 && lat < 53.55);
    }
    /**
     * 经度转换
     * @param lon
     * @param lat
     * @return
     */
    private double transformlon(double lon, double lat){
        double ret = 300.0 + lon + 2.0 * lat + 0.1 * lon * lon + 0.1 * lon * lat + 0.1 * Math.sqrt(Math.abs(lon));
        ret += (20.0 * Math.sin(6.0 * lon * PI) + 20.0 * Math.sin(2.0 * lon * PI)) * 2.0 / 3.0;
        ret += (20.0 * Math.sin(lon * PI) + 40.0 * Math.sin(lon / 3.0 * PI)) * 2.0 / 3.0;
        ret += (150.0 * Math.sin(lon / 12.0 * PI) + 300.0 * Math.sin(lon / 30.0 * PI)) * 2.0 / 3.0;
        return ret;
    }
    /**
     * 纬度转换
     * @param lon
     * @param lat
     * @return
     */
    private double transformlat(double lon, double lat){
        double ret = -100.0 + 2.0 * lon + 3.0 * lat + 0.2 * lat * lat + 0.1 * lon * lat + 0.2 * Math.sqrt(Math.abs(lon));
        ret += (20.0 * Math.sin(6.0 * lon * PI) + 20.0 * Math.sin(2.0 * lon * PI)) * 2.0 / 3.0;
        ret += (20.0 * Math.sin(lat * PI) + 40.0 * Math.sin(lat / 3.0 * PI)) * 2.0 / 3.0;
        ret += (160.0 * Math.sin(lat / 12.0 * PI) + 320 * Math.sin(lat * PI / 30.0)) * 2.0 / 3.0;
        return ret;
    }


    public static void main(String[] args){
        ProjTransform proj = new ProjTransform();
        double[] result = proj.wgs84togcj02(73.6770723800199,39.13769693730567);
        System.out.println(result[0]+", "+result[1]);
    }
}

2.WebProjTrans.java

package com.lzugis.geotools;

import com.lzugis.geotools.utils.ProjTransform;
import com.vividsolutions.jts.geom.*;
import org.geotools.data.FeatureWriter;
import org.geotools.data.FileDataStoreFactorySpi;
import org.geotools.data.Transaction;
import org.geotools.data.shapefile.ShapefileDataStore;
import org.geotools.data.shapefile.ShapefileDataStoreFactory;
import org.geotools.data.simple.SimpleFeatureIterator;
import org.geotools.data.simple.SimpleFeatureSource;
import org.geotools.feature.simple.SimpleFeatureTypeBuilder;
import org.geotools.geometry.jts.JTSFactoryFinder;
import org.geotools.referencing.CRS;
import org.opengis.feature.simple.SimpleFeature;
import org.opengis.feature.simple.SimpleFeatureType;
import org.opengis.referencing.crs.CoordinateReferenceSystem;

import java.io.File;
import java.io.Serializable;
import java.util.HashMap;
import java.util.Map;

/**
 * 实现wgs84，bd09，gcj02三者坐标间的转换
 */
public class WebProjTrans {
    private ProjTransform proj = new ProjTransform();

    public void transformShp(String inputShp, String outputShp, String inCrs, String outCrs){
        try {
            ShapefileDataStore shapeDS = (ShapefileDataStore) new ShapefileDataStoreFactory().createDataStore(new File(inputShp).toURI().toURL());
            SimpleFeatureSource fs = shapeDS.getFeatureSource(shapeDS.getTypeNames()[0]);
            SimpleFeatureIterator it = fs.getFeatures().features();

            FileDataStoreFactorySpi factory = new ShapefileDataStoreFactory();

            Map<String, Serializable> params = new HashMap<String, Serializable>();
            params.put(ShapefileDataStoreFactory.URLP.key, new File(outputShp).toURI().toURL());
            ShapefileDataStore ds = (ShapefileDataStore) factory.createNewDataStore(params);

            CoordinateReferenceSystem crs = CRS.decode("EPSG:4326");
            ds.createSchema(SimpleFeatureTypeBuilder.retype(fs.getSchema(), crs));

            FeatureWriter<SimpleFeatureType, SimpleFeature> writer = ds.getFeatureWriter(ds.getTypeNames()[0], Transaction.AUTO_COMMIT);

            while (it.hasNext()){
                SimpleFeature f = it.next();
                Geometry inGeom = (Geometry)f.getAttribute("the_geom");
                SimpleFeature fNew = writer.next();
                fNew.setAttributes(f.getAttributes());
                Geometry outGeom = projGeometry(inGeom, inCrs, outCrs);
                fNew.setAttribute("the_geom", outGeom);
            }
            writer.write();
            writer.close();
            ds.dispose();
            shapeDS.dispose();
        }catch (Exception e){
            e.printStackTrace();
        }
    }

    /**
     * 点坐标转换
     * @param inCoord
     * @param inCrs
     * @param outCrs
     * @return
     */
    public double[] projPoint(Coordinate inCoord, String inCrs, String outCrs){
        double[] lonlat = new double[2];
        String projStr = inCrs+","+outCrs;
        switch (projStr){
            case "wgs84,bd09":{
                lonlat = proj.wgs84tobd09(inCoord.x, inCoord.y);
                break;
            }
            case "wgs84,gcj02":{
                lonlat = proj.wgs84togcj02(inCoord.x, inCoord.y);
                break;
            }
            case "bd09,wgs84":{
                lonlat = proj.bd09towgs84(inCoord.x, inCoord.y);
                break;
            }
            case "gcj02,wgs84":{
                lonlat = proj.gcj02towgs84(inCoord.x, inCoord.y);
                break;
            }
            case "gcj02,bd09":{
                lonlat = proj.gcj02tobd09(inCoord.x, inCoord.y);
                break;
            }
            case "bd09,gcj02":{
                lonlat = proj.bd09togcj02(inCoord.x, inCoord.y);
                break;
            }
        }
        return lonlat;
    }

    /**
     * 获取转换后的点集合
     * @param geom
     * @param inCrs
     * @param outCrs
     * @return
     */
    public Coordinate[] getProjCoords(Geometry geom, String inCrs, String outCrs){
        Coordinate[] inCoords = geom.getCoordinates();
        Coordinate[] outCoords = new Coordinate[inCoords.length];
        for(int i=0;i<inCoords.length;i++){
            Coordinate inCoord = inCoords[i];
            double[] lonlat = projPoint(inCoord, inCrs, outCrs);
            Coordinate outCoord = new Coordinate(lonlat[0], lonlat[1]);
            String content = lonlat[0]+","+lonlat[1]+";"+inCoord.x+","+inCoord.y;
            outCoords[i] = outCoord;
        }
        return  outCoords;
    }

    /**
     * 单个Geometry做坐标转换
     * @param inGeom
     * @param inCrs
     * @param outCrs
     * @return
     */
    public Geometry projGeometry(Geometry inGeom, String inCrs, String outCrs){
        int geomNum = inGeom.getNumGeometries();
        Geometry[] geoms = new Geometry[geomNum];
        GeometryFactory geometryFactory = JTSFactoryFinder.getGeometryFactory( null );
        for(int i=0;i<geomNum;i++){
            Geometry _inGeom = inGeom.getGeometryN(i),
                    _outGeom = null;
            Coordinate[] _outCoords = getProjCoords(_inGeom, inCrs, outCrs);
            String _geomType = _inGeom.getGeometryType().toLowerCase();
            switch (_geomType){
                case "point":{
                    _outGeom = geometryFactory.createPoint(_outCoords[0]);
                    break;
                }
                case "linestring":{
                    _outGeom = geometryFactory.createLineString(_outCoords);
                    break;
                }
                case "polygon":{
                    _outGeom = geometryFactory.createPolygon(_outCoords);
                    break;
                }
            }
            geoms[i] = _outGeom;
        }
        String geomType = inGeom.getGeometryType().toLowerCase();
        if(geomType.indexOf("multi")>0){//复杂图形
            Geometry outGeom = null;
            if(geomType.indexOf("linestring")>0){//复杂线
                outGeom = geometryFactory.createMultiLineString((LineString[]) geoms);
            }else if(geomType.indexOf("polygon")>0){//复杂面
                outGeom = geometryFactory.createMultiPolygon((Polygon[]) geoms);
            }else{//复杂点
                outGeom = geometryFactory.createMultiPoint((Point[]) geoms);
            }
            return outGeom;
        }else{
            return geoms[0];
        }
    }

    public static void main(String[] args){
        WebProjTrans webProj = new WebProjTrans();
        long start = System.currentTimeMillis();
        String inputShp = "D:\\data\\wgs84\\province.shp",
                outputShp = "D:\\data\\wgs84\\province_gcj02.shp",
                inCrs = "wgs84",//bd09, gcj02
                outCrs = "gcj02";
        webProj.transformShp(inputShp, outputShp, inCrs, outCrs);
        System.out.println("共耗时"+(System.currentTimeMillis() - start)+"ms");
    }
}

---

技术博客
CSDN：http://blog.csdn.NET/gisshixisheng
博客园：http://www.cnblogs.com/lzugis/
在线教程
https://edu.csdn.net/course/detail/7471
Github
https://github.com/lzugis/
联系方式

类型	内容
qq	1004740957
公众号	lzugis15
e-mail	niujp08@qq.com
webgis群	452117357
Android群	337469080
GIS数据可视化群	458292378