GeoTools入门（五）-- CRS操作

最新推荐文章于 2024-06-07 16:22:55 发布

surpassLiang

最新推荐文章于 2024-06-07 16:22:55 发布

阅读量7k

点赞数 12

分类专栏：地理信息 Java 文章标签： geotools CRS EPSG EPSG4326

本文链接：https://blog.csdn.net/oYinHeZhiGuang/article/details/118315724

版权

Java 同时被 2 个专栏收录

60 篇文章 10 订阅

订阅专栏

地理信息

26 篇文章 81 订阅

订阅专栏

1. CRS概述

CRS全称为CoordinateReferenceSystems，中文叫坐标参考系统。它是GeoTools中的一个类，也是我们这次讨论的重点。

在前面我们谈到了JTS库，了解到了数据模型，在这个基础上我们定义了几何体。那么什么是几何体呢？长方形、正方形、长方体等等等等，我们可以求他的周长、面积、体积。然而，在实际生活中，他们仅仅是我们想象的一个物体，并没有任何实际意义。只有当你给他附加一个单位，那么它才具备实际的意义。然而，对于一个地理几何体，我们还需要知道他们的位置。那么记录这些信息的数据结构称之为坐标参考系统。说白了，坐标参考系统就是一个数据结构。在坐标参考系统中，他为我们提供了几个概念：

它定义了使用的轴以及测量单位。因此，您可以从赤道以度为单位测量纬度，从格林威治子午线以度为单位测量经度。或者您可以以米为单位测量 x，以米为单位测量 y，这对于计算距离或面积非常方便。
它定义了世界的形状。不是真的 - 并非所有坐标参考系统都想象世界的相同形状。谷歌使用的 CRS 假装世界是一个完美的球体，而“EPSG:4326”使用的 CRS 有不同的形状——所以如果你把它们混在一起，你的数据将被绘制在错误的地方。这导致的结果就是几何图形不能正确的显示。那是一个很严重的问题。

所以对于任何一个shp文件或者tiff文件，如果他脱离的坐标参考系统或者没有坐标参考系统。他就会变一个毫无意义的几何体或者图片文件。

2. EPSG代码

说到坐标参考系统，我们少不了EPSG，他的全称为欧洲石油标准组织 (EPSG)。这个组织是第一个关心这些东西以将其以数据库形式记录下来的组织。该数据库分布在 Microsoft Access 中，并被移植到各种其他形式。在GeoTools中，它被附带到gt-epsg-hsql的jar中。下面列举一下常见的EPSG代码，如果你想了解更多的EPSG代码，可以访问epsg.io。

2.1 EPSG:4326

EPSG:4326 就是 WGS84 的代码，使用十进制度数通过纬度/经度测量的信息，是目前最流行的地理坐标系统，GPS是基于WGS84的，所以通常我们得到的坐标数据都是WGS84的。一般我们在存储数据时，仍然按WGS84存储。

在这里插入图片描述

2.1.1 WKT定义

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"]]

2.1.2 GeoTools定义坐标

//方式一
CRS.decode("EPSG:4326");
//方式二
DefaultGeographicCRS.WGS84;

GeoTools对于EPSG:4326支持一种默认的定义方式，足见他的重要性了。

2.2 EPSG: 3785

伪墨卡托投影，也被称为球体墨卡托，Web Mercator。它是基于墨卡托投影的，把 WGS84坐标系投影到正方形。我们前面已经知道 WGS84 是基于椭球体的，但是伪墨卡托投影把坐标投影到球体上，这导致两极的失真变大，但是却更容易计算。这也许是为什么被称为”伪“墨卡托吧。另外，伪墨卡托投影还切掉了南北85.051129°纬度以上的地区，以保证整个投影是正方形的。因为墨卡托投影等正形性的特点，在不同层级的图层上物体的形状保持不变，一个正方形可以不断被划分为更多更小的正方形以显示更清晰的细节。很明显，伪墨卡托坐标系是非常显示数据，但是不适合存储数据的，通常我们使用WGS84 存储数据，使用伪墨卡托显示数据。Web Mercator 最早是由 Google 提出的，当前已经成为 Web Map 的事实标准。但是也许是由于上面”伪“的原因，最初 Web Mercator 被拒绝分配EPSG 代码。

在这里插入图片描述

2.2.1 WKT定义

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["X",EAST],
    AXIS["Y",NORTH],
    EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs"],
    AUTHORITY["EPSG","3857"]]

2.2.2 GeoTools定义坐标

CRS.decode("EPSG:3785");

2.3 EPSG:3005

加拿大西海岸的“等面积”投影示例。轴以米为单位测量，便于计算距离或面积。

在这里插入图片描述

2.3.1 WKT定义

PROJCS["NAD83 / BC Albers",
    GEOGCS["NAD83",
        DATUM["North_American_Datum_1983",
            SPHEROID["GRS 1980",6378137,298.257222101,
                AUTHORITY["EPSG","7019"]],
            TOWGS84[0,0,0,0,0,0,0],
            AUTHORITY["EPSG","6269"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4269"]],
    PROJECTION["Albers_Conic_Equal_Area"],
    PARAMETER["standard_parallel_1",50],
    PARAMETER["standard_parallel_2",58.5],
    PARAMETER["latitude_of_center",45],
    PARAMETER["longitude_of_center",-126],
    PARAMETER["false_easting",1000000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH],
    AUTHORITY["EPSG","3005"]]

2.3.2 GeoTools定义坐标

CRS.decode("EPSG:3005");

2.4 EPSG:4490

EPSG:4490为大地坐标系_国家2000大地坐标系CGCS2000。

在这里插入图片描述

2.4.1 WKT定义

GEOGCS["China Geodetic Coordinate System 2000",
    DATUM["China_2000",
        SPHEROID["CGCS2000",6378137,298.257222101,
            AUTHORITY["EPSG","1024"]],
        AUTHORITY["EPSG","1043"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433,
        AUTHORITY["EPSG","9122"]],
    AUTHORITY["EPSG","4490"]]

2.4.2 GeoTools定义坐标

CRS.decode("EPSG:4490");

2.5 EPSG:4479

EPSG:4490为中国大地坐标系2000，它和EPSG:4490的区别在于单位不同，EPSG:4490以度为单位；EPSG:4479以米为单位。

2.5.1 WKT定义

GEOCCS["China Geodetic Coordinate System 2000",
    DATUM["China_2000",
        SPHEROID["CGCS2000",6378137,298.257222101,
            AUTHORITY["EPSG","1024"]],
        AUTHORITY["EPSG","1043"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Geocentric X",OTHER],
    AXIS["Geocentric Y",OTHER],
    AXIS["Geocentric Z",NORTH],
    AUTHORITY["EPSG","4479"]]

2.5.2 GeoTools定义坐标

CRS.decode("EPSG:4479");

3. GeoTools对CRS的操作

GeoTools对CRS的操作主要包括创建、获取和转换。

3.1 CRS的创建

GeoTools对于CRS的创建主要用于创建shp文件。前面的文章里我们有一个示例，通过CSV转换为shp文件，在构造FeatureType的时候需要CRS的创建。而当时我们直接使用DefaultGeographicCRS.WGS84进行创建。当然，上面也提到了也可以通过CRS.decode(“EPSG:4326”); 创建。

3.2 CRS的获取

GeoTools对于CRS的获取主要在解析shp文件的时候用到。我们知道他是创建于FeatureType,在ShapefileDataStore和Geometry都有定义。因此，CRS既可以通过ShapefileDataStore获取，也可以通过Geometry获取。按照正常情况，这两个地方的坐标系统是一致的。

可以通过ShapefileDataStore获取

FileDataStore store = FileDataStoreFinder.getDataStore(sourceFile);
SimpleFeatureType schema = store.getSchema();
CoordinateReferenceSystem coordinateReferenceSystem = schema.getCoordinateReferenceSystem();

通过Geometry获取

//方式一 
Point point = geometryFactory.createPoint(new Coordinate(longitude, latitude));
int srid = point.getSRID();
CoordinateReferenceSystem decode = CRS.decode("EPSG:" + srid);
//方式二
SimpleFeature feature = featureBuilder.buildFeature(null);
feature.getType().getCoordinateReferenceSystem()

3.3 CRS的转换

在实际的应用中，我们还有一种场景，就是不同坐标系统的转换，例如EPSG:4326与EPSG:3857的转换等等。具体代码如下:

CoordinateReferenceSystem dataCRS = schema.getCoordinateReferenceSystem();
CoordinateReferenceSystem worldCRS = map.getCoordinateReferenceSystem();
// allow for some error due to different datums
boolean lenient = true; 
MathTransform transform = CRS.findMathTransform(dataCRS, worldCRS, lenient);
Geometry newGeometry = JTS.transform(oldGeometry, transform);

4. GeoTools对CRS操作示例

通过显示 shapefile 对坐标参考系统进行了直观演示，并展示了更改地图投影如何改变要素的几何形状。

4.1 pom依赖

我们上面提到，EPSG被附带到gt-epsg-hsql的jar中,所以我们需要增加对于gt-epsg-hsql的依赖，具体的pom文件如下：

 <dependencies>
    <dependency>
        <groupId>org.geotools</groupId>
        <artifactId>gt-shapefile</artifactId>
        <version>${geotools.version}</version>
    </dependency>
    <dependency>
        <groupId>org.geotools</groupId>
        <artifactId>gt-swing</artifactId>
        <version>${geotools.version}</version>
    </dependency>
    <dependency>
        <groupId>org.geotools</groupId>
        <artifactId>gt-epsg-hsql</artifactId>
        <version>${geotools.version}</version>
    </dependency>
</dependencies>
<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>

4.2 创建主类CRSLab.java

package cn.surpass.geotools.tutorial.crs;

import java.awt.event.ActionEvent;
import java.io.File;
import java.io.Serializable;
import java.util.HashMap;
import java.util.Map;
import javax.swing.Action;
import javax.swing.JButton;
import javax.swing.JOptionPane;
import javax.swing.JToolBar;
import javax.swing.SwingWorker;
import org.geotools.data.DataStore;
import org.geotools.data.DataStoreFactorySpi;
import org.geotools.data.DefaultTransaction;
import org.geotools.data.FeatureWriter;
import org.geotools.data.FileDataStore;
import org.geotools.data.FileDataStoreFinder;
import org.geotools.data.Query;
import org.geotools.data.Transaction;
import org.geotools.data.shapefile.ShapefileDataStoreFactory;
import org.geotools.data.simple.SimpleFeatureCollection;
import org.geotools.data.simple.SimpleFeatureIterator;
import org.geotools.data.simple.SimpleFeatureSource;
import org.geotools.data.simple.SimpleFeatureStore;
import org.geotools.feature.simple.SimpleFeatureTypeBuilder;
import org.geotools.geometry.jts.JTS;
import org.geotools.map.FeatureLayer;
import org.geotools.map.Layer;
import org.geotools.map.MapContent;
import org.geotools.referencing.CRS;
import org.geotools.styling.SLD;
import org.geotools.styling.Style;
import org.geotools.swing.JMapFrame;
import org.geotools.swing.JProgressWindow;
import org.geotools.swing.action.SafeAction;
import org.geotools.swing.data.JFileDataStoreChooser;
import org.locationtech.jts.geom.Geometry;
import org.opengis.feature.Feature;
import org.opengis.feature.FeatureVisitor;
import org.opengis.feature.simple.SimpleFeature;
import org.opengis.feature.simple.SimpleFeatureType;
import org.opengis.feature.type.FeatureType;
import org.opengis.referencing.crs.CoordinateReferenceSystem;
import org.opengis.referencing.operation.MathTransform;
import org.opengis.util.ProgressListener;
/**
 * @author SurpassLiang
 * @version 1.0
 * @ClassName cn.surpass.geotools.tutorial.crs.CRSLab.java
 * @date 2021/6/25 18:42
 * @desc
 *    启动主类
 */
public class CRSLab {
    private File sourceFile;
    private SimpleFeatureSource featureSource;
    private MapContent map;
}

4.3 创建JMapFrame

我们JMapFrame通过向工具栏添加两个按钮来自定义：一个用于检查要素几何是否有效（例如多边形边界已关闭），另一个用于导出重新投影的要素数据。

private void displayShapefile() throws Exception {
    sourceFile = JFileDataStoreChooser.showOpenFile("shp", null);
    if (sourceFile == null) {
        return;
    }
    FileDataStore store = FileDataStoreFinder.getDataStore(sourceFile);
    featureSource = store.getFeatureSource();

    // Create a map context and add our shapefile to it
    map = new MapContent();
    Style style = SLD.createSimpleStyle(featureSource.getSchema());
    Layer layer = new FeatureLayer(featureSource, style);
    map.layers().add(layer);

    // Create a JMapFrame with custom toolbar buttons
    JMapFrame mapFrame = new JMapFrame(map);
    mapFrame.enableToolBar(true);
    mapFrame.enableStatusBar(true);

    JToolBar toolbar = mapFrame.getToolBar();
    toolbar.addSeparator();
    toolbar.add(new JButton(new ValidateGeometryAction()));
    toolbar.add(new JButton(new ExportShapefileAction()));

    // Display the map frame. When it is closed the application will exit
    mapFrame.setSize(800, 600);
    mapFrame.setVisible(true);
}

4.4 验证几何

我们的工具栏操作是作为嵌套类实现的，大部分工作由父类中的辅助方法完成。

创建ValidateGeometryAction上一节中提到的作为内部类。

class ValidateGeometryAction extends SafeAction {
  ValidateGeometryAction() {
      super("Validate geometry");
      putValue(Action.SHORT_DESCRIPTION, "Check each geometry");
  }

  @Override
  public void action(ActionEvent e) throws Throwable {
      int numInvalid = validateFeatureGeometry(null);
      String msg;
      if (numInvalid == 0) {
          msg = "All feature geometries are valid";
      } else {
          msg = "Invalid geometries: " + numInvalid;
      }
      JOptionPane.showMessageDialog(
          null, msg, "Geometry results", JOptionPane.INFORMATION_MESSAGE);
  }
}

此方法会检查与 shapefile 中每个要素相关联的几何图形是否存在常见问题（例如没有闭合边界的多边形），下面的代码放到上面的内部类中。

private int validateFeatureGeometry(ProgressListener progress) throws Exception {
  final SimpleFeatureCollection featureCollection = featureSource.getFeatures();

  // Rather than use an iterator, create a FeatureVisitor to check each fature
  class ValidationVisitor implements FeatureVisitor {
      public int numInvalidGeometries = 0;

      @Override
      public void visit(Feature f) {
          SimpleFeature feature = (SimpleFeature) f;
          Geometry geom = (Geometry) feature.getDefaultGeometry();
          if (geom != null && !geom.isValid()) {
              numInvalidGeometries++;
              System.out.println("Invalid Geoemtry: " + feature.getID());
          }
      }
  }

  ValidationVisitor visitor = new ValidationVisitor();

  // Pass visitor and the progress bar to feature collection
  featureCollection.accepts(visitor, progress);
  return visitor.numInvalidGeometries;
}

4.5 导出重新投影的形状文件

下一步将形成一个小实用程序，它可以读取 shapefile 并在不同的坐标参考系统中写出 shapefile。从这个实验室中学到的一件重要事情是 MathTransform在两个CoordinateReferenceSystems. 您可以使用 MathTransform一次转换一个点；或使用JTS实用程序类创建一个Geometry修改点的副本。我们使用类似的步骤来导出 CSV2SHAPE 示例中使用的 shapefile。在本例中，我们使用FeatureIterator;从现有 shapefile 中读取内容。并使用FeatureWriter.一次写出一个内容。使用后请关闭这些物品。

action 是一个嵌套类，它委托给exportToShapefile 父类中的方法。

class ExportShapefileAction extends SafeAction {
  ExportShapefileAction() {
      super("Export...");
      putValue(Action.SHORT_DESCRIPTION, "Export using current crs");
  }

  public void action(ActionEvent e) throws Throwable {
      exportToShapefile();
  }
}

将重新投影的数据导出到 shapefile：

private void exportToShapefile() throws Exception {
  SimpleFeatureType schema = featureSource.getSchema();
  JFileDataStoreChooser chooser = new JFileDataStoreChooser("shp");
  chooser.setDialogTitle("Save reprojected shapefile");
  chooser.setSaveFile(sourceFile);
  int returnVal = chooser.showSaveDialog(null);
  if (returnVal != JFileDataStoreChooser.APPROVE_OPTION) {
      return;
  }
  File file = chooser.getSelectedFile();
  if (file.equals(sourceFile)) {
      JOptionPane.showMessageDialog(null, "Cannot replace " + file);
      return;
  }

设置用于处理数据的数学变换：

CoordinateReferenceSystem dataCRS = schema.getCoordinateReferenceSystem();
CoordinateReferenceSystem worldCRS = map.getCoordinateReferenceSystem();
boolean lenient = true; // allow for some error due to different datums
MathTransform transform = CRS.findMathTransform(dataCRS, worldCRS, lenient);

获取所有要素：

SimpleFeatureCollection featureCollection = featureSource.getFeatures();

要创建一个新的 shapefile，我们需要生成一个FeatureType与我们原来的类似的文件。唯一的区别是 CoordinateReferenceSystem几何描述符的。

DataStoreFactorySpi factory = new ShapefileDataStoreFactory();
Map<String, Serializable> create = new HashMap<>();
create.put("url", file.toURI().toURL());
create.put("create spatial index", Boolean.TRUE);
DataStore dataStore = factory.createNewDataStore(create);
SimpleFeatureType featureType = SimpleFeatureTypeBuilder.retype(schema, worldCRS);
dataStore.createSchema(featureType);
// Get the name of the new Shapefile, which will be used to open the FeatureWriter
String createdName = dataStore.getTypeNames()[0];

我们现在可以小心地打开一个迭代器来遍历内容，并打开一个编写器来写出新的 Shapefile。

Transaction transaction = new DefaultTransaction("Reproject");
try (FeatureWriter<SimpleFeatureType, SimpleFeature> writer =
   dataStore.getFeatureWriterAppend(createdName, transaction);
   SimpleFeatureIterator iterator = featureCollection.features()) {
  while (iterator.hasNext()) {
      // copy the contents of each feature and transform the geometry
      SimpleFeature feature = iterator.next();
      SimpleFeature copy = writer.next();
      copy.setAttributes(feature.getAttributes());

      Geometry geometry = (Geometry) feature.getDefaultGeometry();
      Geometry geometry2 = JTS.transform(geometry, transform);

      copy.setDefaultGeometry(geometry2);
      writer.write();
  }
  transaction.commit();
  JOptionPane.showMessageDialog(null, "Export to shapefile complete");
} catch (Exception problem) {
  problem.printStackTrace();
  transaction.rollback();
  JOptionPane.showMessageDialog(null, "Export to shapefile failed");
} finally {
  transaction.close();
}

4.6 运行应用程序

4.6.1 要在地图投影之间切换：

启动应用程序时，系统会提示您输入要显示的 shapefile。在下面的屏幕截图中，我们使用了bc_border地图。
GeoTools 包括一个非常广泛的由 EPSG 参考编号定义的地图投影数据库。对于我们的示例 shapefile，一个合适的替代地图投影是BC Albers。您可以通过键入3005在选择器列表中快速找到它。
请注意，当您将鼠标移到地图上时，坐标现在以米（适用于BC Albers 投影的测量单位）而不是度数显示。
要返回原始投影，请再次打开 CRS 选择器并键入，默认地理投影为4326。

4.6.2 导出重新投影的数据：

当您更改显示的地图投影时，shapefile 中的数据保持不变。

使用bc_border形状文件，特征数据仍然以度为单位，但是当我们选择BC Albers投影时，特征会即时重新投影。
设置重投影数据的显示（例如bc_border shapefile 的3005 BC Albers ）。
单击验证几何按钮以检查要素几何是否正常。
如果没有几何问题，请单击“*导出”*按钮并输入新 shapefile 的名称和路径。

surpassLiang

关注

12
点赞
踩
23

收藏

觉得还不错? 一键收藏
打赏
4
评论
GeoTools入门（五）-- CRS操作

1. CRS概述CRS全称为CoordinateReferenceSystems，中文叫坐标参考系统。它是GeoTools中的一个类，也是我们这次讨论的重点。在前面我们谈到了JTS库，了解到了数据模型，在这个基础上我们定义了几何体。那么什么是几何体呢？长方形、正方形、长方体等等等等，我们可以求他的周长、面积、体积。然而，在实际生活中，他们仅仅是我们想象的一个物体，并没有任何实际意义。只有当你给他附加一个单位，那么它才具备实际的意义。然而，对于一个地理几何体，我们还需要知道他们的位置。那么记录这些信息的数
复制链接

扫一扫