Java geotools 坐标转换

点坐标转换

我点坐标转换为4326到3857，虽然我平时都是4546到3857，都一样。主要时因为我手里没有4546的坐标，

pom文件

这个geotools的依赖比较“皮干”，要是下载不成功，就去官网geotools下的quickstart 的Mavenquick quickstart下面找仓库地址。

  <repositories><repository><id>osgeo</id><name>OSGeo Release Repository</name><url>https://repo.osgeo.org/repository/release/</url><snapshots><enabled>false</enabled></snapshots><releases><enabled>true</enabled></releases></repository><repository><id>osgeo-snapshot</id><name>OSGeo Snapshot Repository</name><url>https://repo.osgeo.org/repository/snapshot/</url><snapshots><enabled>true</enabled></snapshots><releases><enabled>false</enabled></releases></repository></repositories><dependency><groupId>org.geotools</groupId><artifactId>gt-epsg-hsql</artifactId><version>20.0</version></dependency>
```<dependency><groupId>org.geotools</groupId><artifactId>gt-referencing</artifactId><version>20.0</version></dependency><dependency><groupId>org.geotools</groupId><artifactId>gt-api</artifactId><version>20.0</version></dependency>

代码

package com.wangkang.gts4vect;import org.geotools.geometry.jts.JTS;import org.geotools.referencing.CRS;
import org.locationtech.jts.geom.Coordinate;import org.opengis.referencing.FactoryException;
import org.opengis.referencing.crs.CoordinateReferenceSystem;
import org.opengis.referencing.operation.MathTransform;import java.text.DecimalFormat;public class MyGeotools {public static void main(String[] args) {//        GeometryFactory geometryFactory = JTSFactoryFinder.getGeometryFactory();
//        WKTReader reader = new WKTReader( geometryFactory );//某点的经纬度Double y = 34.26095856;Double x = 108.94237668;Double[] coordinate = getCoordinate(x, y);System.out.println(coordinate[0] + " " + coordinate[1]);}//坐标转换public static Double[] getCoordinate(Double x, Double y) {Double[] res = new Double[2];Coordinate tar = null;try {//封装点，这个是通用的，也可以用POINT（y,x）// private static WKTReader reader = new WKTReader( geometryFactory );Coordinate sour = new Coordinate(y, x);//这里要选择转换的坐标系是可以随意更换的CoordinateReferenceSystem source = CRS.decode("EPSG:4326");CoordinateReferenceSystem target = CRS.decode("EPSG:3857");//建立转换，下面两个我屏掉的转换方式会报出需要3/7参数的异常// MathTransform mathTransform = CRS.findMathTransform(source, target);//MathTransform mathTransform1 = CRS.findMathTransform(source, target, false);MathTransform transform = CRS.findMathTransform(source, target, true);tar = new Coordinate();//转换JTS.transform(sour, tar, transform);} catch (FactoryException | org.opengis.referencing.operation.TransformException e) {e.printStackTrace();}String[] split = (tar.toString().substring(1, tar.toString().length() - 1)).split(",");//经纬度精度DecimalFormat fm = new DecimalFormat("0.0000000");res[0] = Double.valueOf(fm.format(Double.valueOf(split[0])));res[1] = Double.valueOf(fm.format(Double.valueOf(split[1])));return res;}
}

结果如下

在这个转换类中，x,y 的位置调换了，这个是我们国家和国外的差异，约定成俗的东西，记住就好，为什么呢？因为不记住接下来的shp坐标系转换图像就飘了，飘到。。。。
如果说整个科学技术不好看，好吧，我再换种写法

package com.wangkang.gts4vect;import org.geotools.geometry.jts.JTS;
import org.geotools.geometry.jts.JTSFactoryFinder;
import org.geotools.referencing.CRS;
import org.locationtech.jts.geom.Coordinate;
import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.geom.GeometryFactory;
import org.locationtech.jts.io.ParseException;
import org.locationtech.jts.io.WKTReader;
import org.opengis.referencing.FactoryException;
import org.opengis.referencing.crs.CoordinateReferenceSystem;
import org.opengis.referencing.operation.MathTransform;
import org.opengis.referencing.operation.TransformException;import java.text.DecimalFormat;public class MyGeotools2 {static  GeometryFactory geometryFactory = JTSFactoryFinder.getGeometryFactory();static  WKTReader reader = new WKTReader( geometryFactory );public static void main(String[] args) {//某点的经纬度Double y=34.26095856;Double x=108.94237668;try {//封装点，x,y 要反过来，不反过来就报纬度超过90的异常String point="POINT("+y+" "+x+")";Geometry geometry = reader.read(point);//这里要选择转换的坐标系是可以随意更换的CoordinateReferenceSystem source = CRS.decode("EPSG:4326");CoordinateReferenceSystem target = CRS.decode("EPSG:3857");//建立转换，下面两个我屏掉的转换方式会报出需要3/7参数的异常// MathTransform mathTransform = CRS.findMathTransform(source, target);//MathTransform mathTransform1 = CRS.findMathTransform(source, target, false);MathTransform transform = CRS.findMathTransform(source, target,true);//转换Geometry transform1 = JTS.transform(geometry, transform);System.out.println(transform1);} catch (FactoryException | TransformException | ParseException e) {e.printStackTrace();}}}

结果如下

其实整个转化就两句代码
MathTransform transform = CRS.findMathTransform(source, target,true);
Geometry transform1 = JTS.transform(geometry, transform);
就这两句。

shp坐标系转换

依赖自己整，上工具类

package com.wangkang.gts4vect;import org.geotools.data.FeatureWriter;
import org.geotools.data.FileDataStoreFactorySpi;
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.JTS;
import org.geotools.geometry.jts.JTSFactoryFinder;
import org.geotools.referencing.CRS;
import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.geom.GeometryFactory;
import org.locationtech.jts.io.WKTReader;
import org.opengis.feature.simple.SimpleFeature;
import org.opengis.feature.simple.SimpleFeatureType;
import org.opengis.referencing.crs.CoordinateReferenceSystem;
import org.opengis.referencing.operation.MathTransform;
import static org.geotools.data.Transaction.AUTO_COMMIT;import java.io.File;
import java.io.Serializable;
import java.util.HashMap;
import java.util.Map;public class geotools3 {private static GeometryFactory geometryFactory = JTSFactoryFinder.getGeometryFactory( null );private static WKTReader reader = new WKTReader( geometryFactory );/*** 单个geom转换* geotools的安装包中的x,y坐标是反着的，因此，直接用geomtry转换，导出的数据是错误，先的将* x,y的坐标换过来。所以，我选择拼装geomtry。* @param geom* @return*/public static Geometry lonlat2WebMactor(Geometry geom,String epsg){try{//CoordinateReferenceSystem的生成方法有两种，EPSG和WKT，然而，我发现把我需要的坐标系的prj文件中的内容用来转换会报“不存在。。。。什么这样的。。”//所以采用EPSG，即使使用了EPSG，导出的转换后的投影文件的prj仍然不对，arcgis 无法识别。当然，同一个坐标系下的prj 文件的内容是相同的，用别的替换即可。CoordinateReferenceSystem crsresource = CRS.decode("EPSG:"+epsg);CoordinateReferenceSystem crsTarget = CRS.decode("EPSG:3857");//CoordinateReferenceSystem crsTarget = CRS.decode(strWKTMercator);// 投影转换MathTransform transform = CRS.findMathTransform(crsresource, crsTarget,true);String geom1 = geom.toString();String[] split = geom1.substring(geom1.lastIndexOf("(")+1, geom1.length() - 3).split(",");StringBuffer res=new StringBuffer();res.append("MULTIPOLYGON (((");for(int i=0;i<split.length;i++){int space = split[i].lastIndexOf(" ");Double x=Double.valueOf(split[i].substring(0,space).trim());Double y=Double.valueOf(split[i].substring(space).trim());//使用coordinate时，产生的结果dest中的x坐标的位数是科学计数
//                Coordinate source = new Coordinate(y, x);
//                Coordinate dest = new Coordinate();
//                JTS.transform(source,dest, transform);
//                String[] str = dest.toString().substring(1, dest.toString().length() - 1).split(",");String point="POINT("+y+" "+x+")";Geometry pointGeom = reader.read(point);Geometry resGeom = JTS.transform(pointGeom, transform);String[] str = resGeom.toString().substring(resGeom.toString().lastIndexOf("(")+1, resGeom.toString().length() - 1).split("\\s+");res.append(str[0]+" "+str[1]+",");}String result = res.substring(0, res.length() - 1);result=result+")))";Geometry read = reader.read(result);return read;}catch (Exception e) {e.printStackTrace();return null;}}/*** 坐标转换：WGS84转墨卡托* shp文件坐标转换。本质是每一个geom转换* @param inputShp shp输入路径* @param outputShp shp输出路径     * * @param  epsg 输入文件的epsg*/public void projectShape(String inputShp, String outputShp,String epsg){try {//源shape文件ShapefileDataStore shapeDS = (ShapefileDataStore) new ShapefileDataStoreFactory().createDataStore(new File(inputShp).toURI().toURL());//创建目标shape文件对象Map<String, Serializable> params = new HashMap<String, Serializable>();FileDataStoreFactorySpi factory = new ShapefileDataStoreFactory();params.put(ShapefileDataStoreFactory.URLP.key, new File(outputShp).toURI().toURL());ShapefileDataStore ds = (ShapefileDataStore) factory.createNewDataStore(params);// 设置属性SimpleFeatureSource fs = shapeDS.getFeatureSource(shapeDS.getTypeNames()[0]);//CoordinateReferenceSystem crs = CRS.parseWKT(strWKTMercator);CoordinateReferenceSystem crs = CRS.decode("EPSG:3857");ds.createSchema(SimpleFeatureTypeBuilder.retype(fs.getSchema(),crs));//设置writerFeatureWriter<SimpleFeatureType, SimpleFeature> writer = ds.getFeatureWriter(ds.getTypeNames()[0],AUTO_COMMIT);//写记录SimpleFeatureIterator it = fs.getFeatures().features();try {while (it.hasNext()) {SimpleFeature f = it.next();SimpleFeature fNew = writer.next();fNew.setAttributes(f.getAttributes());//每个geom转换Geometry geom = lonlat2WebMactor((Geometry)f.getAttribute("the_geom"),epsg);fNew.setAttribute("the_geom", geom);}}finally {it.close();}writer.write();writer.close();ds.dispose();shapeDS.dispose();}catch (Exception e) {e.printStackTrace();}}public static void main(String[] args) {}}

上面的只能转换面，写的也很蹩脚，令人痛心疾首，所以lonlat2WebMactor方法更新如下

   try{        CoordinateReferenceSystem crsresource = CRS.decode("EPSG:"+epsg);CoordinateReferenceSystem crsTarget = CRS.decode("EPSG:3857");// 投影转换MathTransform transform = CRS.findMathTransform(crsresource, crsTarget,true);Coordinate[] coordinates = geom.getCoordinates();for (Coordinate coordinate : coordinates) {               double x = coordinate.getX();coordinate.setX(coordinate.getY());coordinate.setY(x);}Geometry transform1 = JTS.transform(geom, transform);return transform1;}catch (Exception e) {e.printStackTrace();return null;}

转换的功能已经完成了。
需要注意的是：
1.转换完成后(我这里是3857)，prj的文件内容是有问题的。arcgis，QGIS并不能识别这个空间坐标系
如图：有问题的

用arcgis转换出的3857的prj是：

String prj="PROJCS[\"WGS_1984_Web_Mercator_Auxiliary_Sphere\","+ "GEOGCS[\"GCS_WGS_1984\","+ "DATUM[\"D_WGS_1984\","+ "SPHEROID[\"WGS_1984\",6378137.0,298.257223563]],"+ "PRIMEM[\"Greenwich\",0.0],"+ "UNIT[\"Degree\",0.0174532925199433]],"+ "PROJECTION[\"Mercator_Auxiliary_Sphere\"],"+ "PARAMETER[\"False_Easting\",0.0],"+ "PARAMETER[\"False_Northing\",0.0],"+ "PARAMETER[\"Central_Meridian\",0.0],"+ "PARAMETER[\"Standard_Parallel_1\",0.0],"+ "PARAMETER[\"Auxiliary_Sphere_Type\",0.0],"+ "UNIT[\"Meter\",1.0]]";

测试阶段，可以手动把prj的内容替换了。我用的时候写了一个扣JIO的东西。目的，替换prj中的内容

package com.suntoon.util;import com.suntoon.mapper.CoordinateSystemMapper;
import org.springframework.beans.factory.annotation.Autowired;
import org.springframework.stereotype.Component;import javax.annotation.Resource;
import java.io.*;
@Component
public class PrjUtil {/*** 修改prj文件中的坐标系* path:prj文件的全路径*/public void editPrj(String path){final String MKT84 = "PROJCS[\"WGS_1984_Web_Mercator_Auxiliary_Sphere\","+ "GEOGCS[\"GCS_WGS_1984\","+ "DATUM[\"D_WGS_1984\","+ "SPHEROID[\"WGS_1984\",6378137.0,298.257223563]],"+ "PRIMEM[\"Greenwich\",0.0],"+ "UNIT[\"Degree\",0.0174532925199433]],"+ "PROJECTION[\"Mercator_Auxiliary_Sphere\"],"+ "PARAMETER[\"False_Easting\",0.0],"+ "PARAMETER[\"False_Northing\",0.0],"+ "PARAMETER[\"Central_Meridian\",0.0],"+ "PARAMETER[\"Standard_Parallel_1\",0.0],"+ "PARAMETER[\"Auxiliary_Sphere_Type\",0.0],"+ "UNIT[\"Meter\",1.0]]";File file = new File(path);try (FileOutputStream outputStream = new FileOutputStream(file);OutputStreamWriter outputStreamWriter = new OutputStreamWriter(outputStream,"utf-8");BufferedWriter bufferedWriter = new BufferedWriter(outputStreamWriter);){bufferedWriter.write(MKT84);} catch (IOException e) {e.printStackTrace();}}}

2.切记x,y是位置换了的，不然，有时候不会报错，改完prj的文件arcgis可以打开，但是，和用arcgis转出的图层的样子，属性内容是一样的，但是坐标不对（可以发布看下返回链接的bbox的经纬度，也可以几何计算看看）。
3.我把用arcgis转换和用geotools转换的结果做了对比
geotools:

arcgis:

但从边界看起来，似乎是有误差，而且误差不小。但是，坐标是一摸一样的。在显示的位数中，没有误差。后来放大图层，观看边界，发现图形是完全重合的。这个显示的问题并不会影响获取属性数据。