GeoTools 自定義坐標系

前言

在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. 開發環境

本文使用如下開發環境，以供參考。
:::block-1
時間：2025年

GeoTools：v34-SNAPSHOT

IDE：IDEA2025.1.2

JDK：v17

OpenLayers：v9.2.4

Layui：v2.9.14
:::

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 formatfinal 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 formatfinal 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 formatfinal 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_34return 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權威機構為例，輸出信息如下，首次調用輸出情況可能會稍慢，數據量還是有點多。
:::block-1
authorityCodes：[AUTO:42001, AUTO:42002, AUTO:42003, AUTO:42004, AUTO:97001,······]

默認支持的EPSG坐標系數量：7921

默認支持的EPSG坐標：[WGS84(DD), 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007,······]
:::
檢查目標坐標系是否存在，如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地理坐標系，元數據定義如下：
:::block-1
GEOGCS["Geographic Coordinate System",DATUM["D_WGS84",SPHEROID["WGS84",6378137,298.257223560493]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]]
:::
（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視圖
:::block-1
SELECT f_table_name, f_geometry_column, srid FROM geometry_columns WHERE f_table_name = ‘你的表名’; // 如 province
:::
該方式也會返回幾何列的空間參考ID（SRID）。