❝

在GIS开发中，坐标系统是重中之重，在接到任务时首先要确定的就是坐标系。大多数地图库或者互联网地图默认支持WGS84地理坐标系和Web墨卡托投影坐标系。而在我国要求使用自然资源数据使用2000国家大地坐标（CGCS2000）。

1. 背景

经国务院批准，我国自2008年7月1日起，启用2000国家大地坐标系。

2000国家大地坐标系与原有国家大地坐标系转换、衔接的过渡期为8至10年。根据原国家测绘地理信息局下发的《关于加快2000国家大地坐标系推广使用的通知》（国测国发﹝2013﹞11号）和原国土资源部、原国家测绘地理信息局联合发文《关于加快使用2000国家大地坐标系的通知》（国土资发﹝2017﹞30号）的文件要求，明确自2018年7月1日起全面使用2000国家大地坐标系。

2018年6月底前完成全系统各类国土资源空间数据向2000国家大地坐标系转换，7月1日后涉及空间坐标的报部审查和备案项目，全部采用2000国家大地坐标系，不再接受非2000系上报的项目报件。

2018年12月，自然资源部发布《关于停止提供1954年北京坐标系和1980西安坐标系基础测绘成果的公告》（2018年第55号），决定自2019年1月1日起，全面停止向社会提供1954年北京坐标系和1980西安坐标系基础测绘成果。

基于上述内容要求，在生产实践中，从数据生产到数据成果入库，都需要使用2000国家大地坐标系。

2. 开发环境

本文使用如下开发环境，以供参考。

3. 自定义坐标系统

在项目目录utils下新建一个工具类CrsLab，使用final定义几个常用坐标系常量。坐标系的定义有多种标准，包括ogc标准、esri标准、Proj4js以及PostGIS等。在GeoTools中，采用符合行业规范的OGC标准格式。其中OGC格式又分为OGC WKT和OGC WKT2，代码中分别表示为epsg4522wkt和epsg4522wkt2，经过测试验证，目前尚未支持OGC WKT2格式。

最后定义一个坐标系解析方法getCRSByCodeNumber以返回目标坐标系统，该方法接收一个数字字符串参数，即EPSG编码，如WGS84编码为4326、web墨卡托编码为3857（又900913或者3785）、CGCS2000编码为4490。将坐标系常量传递给parseWKT解析为目标坐标系。

package org.geotools.utils;

import org.geotools.api.referencing.crs.CoordinateReferenceSystem;
import org.geotools.referencing.CRS;
import org.slf4j.Logger;
import org.slf4j.LoggerFactory;

public class CrsLab {

    private static final Logger log = LoggerFactory.getLogger(CrsLab.class);

    // ogc wkt format
    final static String epsg4490 = "GEOGCS["China Geodetic Coordinate System 2000",n" +
            "    DATUM["China_2000",n" +
            "        SPHEROID["CGCS2000",6378137,298.257222101,n" +
            "            AUTHORITY["EPSG","1024"]],n" +
            "        AUTHORITY["EPSG","1043"]],n" +
            "    PRIMEM["Greenwich",0,n" +
            "        AUTHORITY["EPSG","8901"]],n" +
            "    UNIT["degree",0.0174532925199433,n" +
            "        AUTHORITY["EPSG","9122"]],n" +
            "    AUTHORITY["EPSG","4490"]]";

    // ogc wkt format
    final static String epsg4522wkt = "PROJCS["CGCS2000 / 3-degree Gauss-Kruger zone 34",n" +
            "    GEOGCS["China Geodetic Coordinate System 2000",n" +
            "        DATUM["China_2000",n" +
            "            SPHEROID["CGCS2000",6378137,298.257222101,n" +
            "                AUTHORITY["EPSG","1024"]],n" +
            "            AUTHORITY["EPSG","1043"]],n" +
            "        PRIMEM["Greenwich",0,n" +
            "            AUTHORITY["EPSG","8901"]],n" +
            "        UNIT["degree",0.0174532925199433,n" +
            "            AUTHORITY["EPSG","9122"]],n" +
            "        AUTHORITY["EPSG","4490"]],n" +
            "    PROJECTION["Transverse_Mercator"],n" +
            "    PARAMETER["latitude_of_origin",0],n" +
            "    PARAMETER["central_meridian",102],n" +
            "    PARAMETER["scale_factor",1],n" +
            "    PARAMETER["false_easting",34500000],n" +
            "    PARAMETER["false_northing",0],n" +
            "    UNIT["metre",1,n" +
            "        AUTHORITY["EPSG","9001"]],n" +
            "    AUTHORITY["EPSG","4522"]]";

    final static String epsg4522wkt2 = "PROJCRS["CGCS2000 / 3-degree Gauss-Kruger zone 34",BASEGEOGCRS["China Geodetic Coordinate System 2000",DATUM["China 2000",ELLIPSOID["CGCS2000",6378137,298.257222101,LENGTHUNIT["metre",1]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],ID["EPSG",4490]],CONVERSION["3-degree Gauss-Kruger zone 34",METHOD["Transverse Mercator",ID["EPSG",9807]],PARAMETER["Latitude of natural origin",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8801]],PARAMETER["Longitude of natural origin",102,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["Scale factor at natural origin",1,SCALEUNIT["unity",1],ID["EPSG",8805]],PARAMETER["False easting",34500000,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",0,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["northing (X)",north,ORDER[1],LENGTHUNIT["metre",1]],AXIS["easting (Y)",east,ORDER[2],LENGTHUNIT["metre",1]],USAGE[SCOPE["Cadastre, engineering survey, topographic mapping (large scale)."],AREA["China - between 100°30'E and 103°30'E."],BBOX[21.13,100.5,42.69,103.5]],ID["EPSG",4522]]";

    // esri wkt format
    final static String epsg4522esriwkt =  "PROJCS["CGCS2000_3_Degree_GK_Zone_34",GEOGCS["GCS_China_Geodetic_Coordinate_System_2000",DATUM["D_China_2000",SPHEROID["CGCS2000",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Gauss_Kruger"],PARAMETER["False_Easting",34500000.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",102.0],PARAMETER["Scale_Factor",1.0],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]]";

    final static String epsg3857 = "PROJCS["WGS 84 / Pseudo-Mercator",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Mercator_1SP"],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs"],AUTHORITY["EPSG","3857"]]";

    public static CoordinateReferenceSystem getCRSByCodeNumber(String epsgCode) throws Exception {
        CoordinateReferenceSystem crs = null;
        if(epsgCode.equals("4490")){
            crs = CRS.parseWKT(epsg4490);
        }else if(epsgCode.equals("4522")){
            crs = CRS.parseWKT(epsg4522esriwkt);
        }else if(epsgCode.equals("3857")){
            crs = CRS.parseWKT(epsg3857);
        }
        System.out.println("坐标系统："+CRS.toSRS(crs));
        // 坐标系统：CGCS2000_3_Degree_GK_Zone_34
        return crs;
    }
}

编写完成以上代码已经可以根据系统要求自定义坐标系了。

以下介绍一些拓展信息。GeoTools中有一个接口CRSAuthorityFactory提供了getSupportedCodes(String authority)方法，用于获取指定权威机构（authority）支持的所有坐标参考系统（CRS）的代码。这里的authority参数指的是定义CRS的权威机构的名称，如“EPSG”。

常见的权威机构包括：

可以通过以下方法获取各权威机构坐标定义数量及编号使用情况。

// 获取 CRS 权威工厂
CRSAuthorityFactory authorityFactory = CRS.getAuthorityFactory(true);

// 指定权威机构并获取支持的代码
String authority = "EPSG"; // 这里可以替换为其他权威机构
Set<String> authorityCodes = authorityFactory.getAuthorityCodes(CoordinateReferenceSystem.class);
System.out.println("authorityCodes："+authorityCodes);

Set<String> supportedCodes = CRS.getSupportedCodes(authority);
// 输出默认支持的坐标系
System.out.println("默认支持的EPSG坐标系数量："+supportedCodes.size());
System.out.println("默认支持的EPSG坐标列表："+supportedCodes);

以EPSG权威机构为例，输出信息如下，首次调用输出情况可能会稍慢，数据量还是有点多。

检查目标坐标系是否存在，如CGCS2000_3_Degree_GK_Zone_34（4522）。

// 检查指定坐标是否被支持
String epsgcode = "4522"; // CGCS2000_3_Degree_GK_Zone_34
if (supportedCodes.contains(epsgcode)) {
    System.out.println(authority + "权威机构支持 EPSG 编码为 " + epsgcode + " 的坐标系");
}else{
    System.out.println(authority + " 权威机构暂不支持 EPSG 编码为 " + epsgcode + " 的坐标系");
}

在EPSG权威机构中进行查找，以上代码将会输出EPSG权威机构支持 EPSG 编码为 4522 的坐标系。

将上述变量epsgcode的值改为4490，通过输出结果EPSG权威机构支持 EPSG 编码为 4490 的坐标系可以看到该坐标系也是受到支持的。

以下代码可以用于获取常用权威机构支持的坐标信息。

// 获取常用权威机构支持的代码
Set<String> epsgCodes = getSupportedCodes("EPSG");  // 标准地理坐标系
Set<String> esriCodes = getSupportedCodes("ESRI");  // Esri特定坐标系
Set<String> iauCodes = getSupportedCodes("IAU");    // 天体坐标系
Set<String> autoCodes = getSupportedCodes("AUTO");  // 自动投影

我们知道在GeoTools中可以直接使用CRS.decode("EPSG:4326")此种方式定义目标坐标系，这种方式使用起来简单明了，那么为什么不直接使用此种方式，而是要自定义坐标系统呢？

下面来看一下我使用两种方式将Shp数据经过坐标转换后导入PostGIS空间数据库的输出结果。坐标系为WGS84地理坐标系，元数据定义如下：

（1）CrsLab.getCRSByCodeNumber("4522")

（2）CRS.decode("EPSG:4522")

从以上输出结果来看，使用两种方式转换的数据，坐标值存在一定的差异。第二种方式x和y值的顺序相反，与工作实际使用不太相符，要怎么解决这种情况呢？需要在坐标定义时多加一个布尔值参数，CoordinateReferenceSystem targetCrs = CRS.decode("EPSG:4522",Boolean.TRUE)，并将参数值设置为TRUE即可。

在GIS软件ArcMap中打开经过投影转换后的Shp数据，其坐标显示结果与使用GeoTools转换入库后的结果保持一致。

4. 坐标转换

在进行坐标转换前，需要获取源数据坐标系信息，然后使用以下任意一种方式定义目标坐标系，使用第二种方式的时候需要注意x和y坐标的顺序，将第二个参数设置为TRUE可使得数据与日常使用习惯保持一致。

坐标转换流程大致如下

// 源数据坐标系
CoordinateReferenceSystem originalCrs = dataSchema.getCoordinateReferenceSystem();
// 方式一：自定义坐标系
// CoordinateReferenceSystem targetCrs =  CrsLab.getCRSByCodeNumber("4522");

// 方式二：使用epsg code 定义坐标系，强制使用经度在前、纬度在后的顺序
CoordinateReferenceSystem targetCrs =  CRS.decode("EPSG:4522",Boolean.TRUE);

使用findMathTransform进行坐标转换。第一个参数为原始坐标系，第二个参数为目标坐标系，第三个参数表示即使没有可用于基准偏移的信息，也应该创建数学变换。默认值为false。

// 定义坐标转换
MathTransform transform = CRS.findMathTransform(originCrs,targetCrs,true);

获取到源几何对象之后，调用JTS.transform方法对其进行坐标转换，最后将转换完成的几何对象写入几何字段。

// 源几何对象
Geometry sourceGeometry = (Geometry) feature.getDefaultGeometry();
// 转换后的几何对象
Geometry crsTransformGeometry = JTS.transform(sourceGeometry, transform);
// 设置参考系数字
crsTransformGeometry.setSRID(4522);
// 写入转换后的几何对象
targetFeature.setAttribute("geom", crsTransformGeometry);

5. 查看数据库坐标系

将成果数据导入PostGIS空间数据库之后，怎么知道使用了正确的坐标系了呢？在PostGIS中，有以下几种方式可以查看数据坐标系。

（1）使用ST_SRID函数

SELECT ST_SRID(geom) AS srid FROM province LIMIT 3;

该方法会返回表中第一个几何列的空间参考ID（SRID）。

（2）查询 geometry_columns视图

该方式也会返回几何列的空间参考ID（SRID）。

6. 坐标拓展

从https://epsg.io网站查询出的坐标系信息有多种格式，包括OGC WKT、OGC WKT2以及ESRI WKT等。正常情况下，可以采用OGC WKT、ESRI WKT两种方式定义坐标系。

（1）OGC WKT格式

（2）OGC WKT2格式

（3）ESRI WKT格式

使用时可以直接选中复制坐标信息，也可以点击【Copy Text】按钮或者【Copy Formatted】按钮进行复制。

7. 参考资料

• EPSG 坐标系网站：https://epsg.io

• GeoTools API：https://docs.geotools.org/latest/javadocs/org/geotools/api/referencing/crs/CoordinateReferenceSystem.html

❝

OpenLayers示例数据下载，请在公众号后台回复：ol数据

全国信息化工程师－GIS 应用水平考试资料，请在公众号后台回复：GIS考试

❝

GIS之路公众号已经接入了智能助手，欢迎大家前来提问。

欢迎访问我的博客网站-长谈GIS：http://shanhaitalk.com