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

经纬度的相互转换"/>

结合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->gcjdouble[] bd09_gcj02 = bd09togcj02(bd_lon, bd_lat);//2、gcj->wgs84return 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->gcjdouble[] wgs84_gcj02 = wgs84togcj02(wgs_lon, wgs_lat);//2、gcj->bd09return 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.05return !(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, gcj02outCrs = "gcj02";webProj.transformShp(inputShp, outputShp, inCrs, outCrs);System.out.println("共耗时"+(System.currentTimeMillis() - start)+"ms");}
}

技术博客
CSDN：.NET/gisshixisheng
博客园：/
在线教程

Github
/
联系方式