geotools学习（三）坐标系

作者: MrSwilder | 来源:发表于2019-10-28 15:15 被阅读0次

坐标系

本教程通过显示一个shapefile对坐标参考系统进行了可视化演示，并展示了更改地图投影如何改变要素的几何形状。
1.请确认pox.xml中包含了下面内容

<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>maven2-repository.dev.java.net</id>
          <name>Java.net repository</name>
          <url>http://download.java.net/maven/2</url>
      </repository>
      <repository>
          <id>osgeo</id>
          <name>Open Source Geospatial Foundation Repository</name>
          <url>http://download.osgeo.org/webdav/geotools/</url>
      </repository>
      <repository>
        <snapshots>
          <enabled>true</enabled>
        </snapshots>
        <id>boundless</id>
        <name>Boundless Maven Repository</name>
        <url>http://repo.boundlessgeo.com/main</url>
      </repository>
  </repositories>

2.创建一个包org.geotools.tutorial.crs和类CRSLab.java，并拷贝下面代码：

package org.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;

/** This is a visual example of changing the coordinate reference system of a feature layer. */
public class CRSLab {

    private File sourceFile;
    private SimpleFeatureSource featureSource;
    private MapContent map;

    public static void main(String[] args) throws Exception {
        CRSLab lab = new CRSLab();
        lab.displayShapefile();
    }

3.注意，我们通过在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);
    }

下面是我们如何配置JMapFrame
1.我们启用了一个状态行;它包含一个按钮，允许改变地图坐标参考系统。
2.我们已经启用了工具栏并向其添加了两个操作(我们将在下一节中定义)。

验证几何

我们的工具栏操作是作为一个嵌套类实现的，大部分工作是由父类中的帮助器方法完成的。
1.创建前一节中提到的作为内部类的ValidateGeometryAction。

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

        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);
        }
    }

2.此方法检查与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;

            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;
    }

导出重投影的shapefile

接下来的操作将形成一个小实用程序，它可以读取一个shapefile并在不同的坐标参考系统中写出一个shapefile。

从这个实验室获得的一个重要信息是，在两个CoordinateReferenceSystems之间创建MathTransform是多么容易。可以使用MathTransform一次转换一个点;或者使用JTS实用工具类创建修改了点的几何图形的副本

我们使用与CSV2SHAPE示例类似的步骤来导出shapefile。在本例中，我们使用featuator从现有的shapefile中读取内容;然后使用功能编辑器一次一个地把内容写出来。请在使用后关闭这些物品。

1.该操作是一个嵌套类，它委托给父类中的exportToShapefile方法。

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

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

2.导出重新投影的数据到一个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;
        }

3.设置一个数学转换用于处理数据

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);

4.抓取所有要素

  SimpleFeatureCollection featureCollection = featureSource.getFeatures();

5.为了创建一个新的shapefile，我们需要生成一个与原始的类似的功能类型。唯一的区别是几何描述符的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];

6.现在，我们可以小心地打开一个迭代器来遍历内容，并使用一个写入器来写出新的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();
        }
    }

运行程序

在地图投影之间切换:
1.当您启动应用程序时，将提示您显示一个shapefile。在下面的屏幕截图中，我们使用了bc_border映射，它可以作为 uDig sample data set的一部分下载

image.png
2.geotools包括一个非常广泛的数据库地图投影由EPSG参考编号定义。对于我们的示例shapefile，一个合适的可选映射投影是BC Albers。
你可以通过输入3005在选择器列表中快速找到它。
当你点击确定，地图会显示在新的投影

image.png
3.请注意，当您将鼠标移到地图上时，坐标现在显示为米(适用于BC Albers投影的测量单位)，而不是度。
4.要返回到原始投影，再次打开CRS选择器并键入4326为默认的地理投影
注意，地图坐标现在再次以度数表示。

导出重新投影的数据:

1.当您更改显示的地图投影时，shapefile中的数据将保持不变。
使用bc_border shapefile，特征数据仍然以度为单位，但是当我们选择BC Albers投影时，特征会动态地重新投影。
2.设置重新投影数据的显示(例如bc_border shapefile的3005 BC Albers)。
3.单击Validate geometry按钮以检查特性几何图形是否正确。
4.如果没有几何图形问题，则单击Export按钮并输入新shapefile的名称和路径。

扩展尝试

对于上面的应用程序，可以尝试以下几件事情：
1.请看屏幕底部显示的EPSG:4326和EPSG:3005中的坐标。你应该能够看到一个是用度来测量的，另一个是用米来测量的。
2.做一个按钮，打印出地图坐标参考系统作为人类可读的wkt格式。这与shapefile的prj辅助文件使用的文本格式相同
3.让用户久等是不礼貌的;SwingWorker类是Java的一部分。
用以下代码替换ValidateGeometryAction:

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

        public void action(ActionEvent e) throws Throwable {
            // Here we use the SwingWorker helper class to run the validation routine in a
            // background thread, otherwise the GUI would wait and the progress bar would not be
            // displayed properly
            SwingWorker worker =
                    new SwingWorker<String, Object>() {
                        protected String doInBackground() throws Exception {
                            // For shapefiles with many features its nice to display a progress bar
                            final JProgressWindow progress = new JProgressWindow(null);
                            progress.setTitle("Validating feature geometry");

                            int numInvalid = validateFeatureGeometry(progress);
                            if (numInvalid == 0) {
                                return "All feature geometries are valid";
                            } else {
                                return "Invalid geometries: " + numInvalid;
                            }
                        }

                        protected void done() {
                            try {
                                Object result = get();
                                JOptionPane.showMessageDialog(
                                        null,
                                        result,
                                        "Geometry results",
                                        JOptionPane.INFORMATION_MESSAGE);
                            } catch (Exception ignore) {
                            }
                        }
                    };
            // This statement runs the validation method in a background thread
            worker.execute();
        }
    }

4.访问JTS网站，查找如何简化几何图形。在将其写出来之前，修改示例以简化几何结构——可以尝试几种技术(推荐使用TopologyPreservingSimplifier和DouglasPeuckerSimplifier类)。

这是一种常见的数据准备形式
1.关于几何图形，特别是你自己做的，有一件事是危险的，那就是它们可能是无效的
修复无效的几何图形有许多技巧。一个简单的起点是使用geometry.buffer(0)。使用这个技巧来构建您自己的shapefile数据清理器。

2.代替手工进行所有几何变换的另一种方法是请求所需投影中的数据。
导出方法的这个版本展示了如何使用查询对象来检索重新投影的特性并将它们写入到一个新的shapefile中，而不是像上面那样“手工”转换特性。

 Query query = new Query(typeName);
        query.setCoordinateSystemReproject(map.getCoordinateReferenceSystem());

        SimpleFeatureCollection featureCollection = featureSource.getFeatures(query);

        // And create a new Shapefile with the results
        DataStoreFactorySpi factory = new ShapefileDataStoreFactory();

        Map<String, Serializable> create = new HashMap<>();
        create.put("url", file.toURI().toURL());
        create.put("create spatial index", Boolean.TRUE);
        DataStore newDataStore = factory.createNewDataStore(create);

        newDataStore.createSchema(featureCollection.getSchema());
        Transaction transaction = new DefaultTransaction("Reproject");
        SimpleFeatureStore featureStore;
        featureStore = (SimpleFeatureStore) newDataStore.getFeatureSource(typeName);
        featureStore.setTransaction(transaction);
        try {
            featureStore.addFeatures(featureCollection);
            transaction.commit();
            JOptionPane.showMessageDialog(
                    null,
                    "Export to shapefile complete",
                    "Export",
                    JOptionPane.INFORMATION_MESSAGE);
        } catch (Exception problem) {
            transaction.rollback();
            problem.printStackTrace();
            JOptionPane.showMessageDialog(
                    null, "Export to shapefile failed", "Export", JOptionPane.ERROR_MESSAGE);
        } finally {
            transaction.close();
        }

Geometry

几何学就是地理信息系统的形状。
通常一个特征只有一个几何形状;属性提供了进一步的描述或度量。有时很难把几何看作是另一个属性。如果考虑同一事物有多个表示形式的情况，这将有所帮助。
我们可以以悉尼市为例：
作为一个单一的地点，即作为城市边界的一个点(这样你就可以知道你什么时候在悉尼)，即一个多边形

Point

下面是一个使用(WKT)格式创建点的示例。

GeometryFactory geometryFactory = JTSFactoryFinder.getGeometryFactory(null);

WKTReader reader = new WKTReader(geometryFactory);
Point point = (Point) reader.read("POINT (1 1)");

您还可以直接使用GeometryFactory手动创建一个点。

GeometryFactory geometryFactory = JTSFactoryFinder.getGeometryFactory(null);
Coordinate coord = new Coordinate(1, 1);
Point point = geometryFactory.createPoint(coord);

Line

下面是一个WKT LineString的例子

GeometryFactory geometryFactory = JTSFactoryFinder.getGeometryFactory(null);

WKTReader reader = new WKTReader(geometryFactory);
LineString line = (LineString) reader.read("LINESTRING(0 2, 2 0, 8 6)");

LineString是片段序列，其方式与java String是字符序列相同。
下面是一个使用GeometryFactory的示例。

GeometryFactory geometryFactory = JTSFactoryFinder.getGeometryFactory(null);

Coordinate[] coords  =
new Coordinate[] {new Coordinate(0, 2), new Coordinate(2, 0), new Coordinate(8, 6) };

LineString line = geometryFactory.createLineString(coordinates);

Polygon

在WKT中通过构造一个外环形成一个多边形，然后是一系列的孔洞。

GeometryFactory geometryFactory = JTSFactoryFinder.getGeometryFactory(null);
WKTReader reader = new WKTReader(geometryFactory);
Polygon polygon = (Polygon) reader.read("POLYGON((20 10, 30 0, 40 10, 30 20, 20 10))");

为什么不使用Java Shape呢

Java Shape实际上是非常有用的，它涵盖了上面提到的思想——但是它非常关注绘图。几何学允许我们处理地理信息系统的“信息”部分——我们可以使用它来创建新的几何学和测试几何学之间的关系。

// Create Geometry using other Geometry
Geometry smoke = fire.buffer(10);
Geometry evacuate = cities.intersection(smoke);

// test important relationships
boolean onFire = me.intersects(fire);
boolean thatIntoYou = me.disjoint(you);
boolean run = you.isWithinDistance(fire, 2);

// relationships actually captured as a fancy
// String called an intersection matrix
//
IntersectionMatrix matrix = he.relate(you);
thatIntoYou = matrix.isDisjoint();

// it really is a fancy string; you can do
// pattern matching to sort out what the geometries
// being compared are up to
boolean disjoint = matrix.matches("FF*FF****");

建议您阅读包含有用定义的JTS的javadocs。

坐标参考系

前面我们讨论了JTS库，它为几何图形提供了数据模型。这是真正的火箭科学的地图制作-一个想法的形状和足够的数学做一些有趣的东西。但有一个问题我们还没有回答——几何是什么意思?

你可能认为我在开玩笑，但问题是严肃的。几何就是一堆数学(数学意义上的一组点)。它们本身没有意义。

一个简单的例子是数字3。数字3本身没有意义。只有当你附加一个“单元”时，数字3才能代表一个概念。3米。3英尺。3年。

为了提供一个有意义的几何图形，我们需要知道这些点的意义。我们需要知道它们的位置，而告诉我们这一点的数据结构被称为坐标参考系统。

坐标系为我们定义了两个概念:

1.它定义了使用的轴和度量单位。
所以你可以用度数从赤道测量纬度，用度数从格林威治子午线测量经度。
或者你可以用米来测量x，用米来测量y，这对于计算距离或面积非常方便。
2.它定义了世界的形状。不，它确实是-不是所有的坐标参考系统想象世界的形状是一样的。谷歌使用的CRS假装世界是一个完美的球体，而“EPSG:4326”使用的CRS在头脑中有一个不同的形状——所以如果你把它们混在一起，你的数据将被画在错误的地方。

作为一名程序员，我认为坐标参考系统的想法是一个巧妙的黑客。每个人都在谈论三维空间中的点，而不是一直记录x、y、z，我们在作弊，只记录两个点(lat、lon)，并使用地球形状的模型来计算z。

EPSG代码

欧洲石油标准组织(EPSG)是第一个非常关心这个问题并将其以数据库形式记录下来的组织。该数据库分布在Microsoft Access中，并被移植到各种其他形式，包括GeoTools中包含的gt-hsql jar
EPSG: 4326
EPSG投影4326 - WGS 84

image.png
这是一个大问题:使用十进制度测量纬度/经度的信息。

CRS.decode("EPSG:4326");

DefaultGeographicCRS.WGS84;

EPSG: 3785
流行的可视化CRS / Mercator

image.png
许多web地图应用程序使用的“谷歌地图”投影的官方代码。假装世界是一个球体是很好的(它使你的数学计算非常快)。然而，如果你在奥斯陆画一个正方形，它看起来真的很奇怪。

CRS.decode("EPSG:3785");

EPSG:3005
NAD83 / BC Albers

image.png
加拿大西海岸“等面积”投影的例子。轴线以米为单位，便于计算距离或面积。

CRS.decode("EPSG:3005");

轴顺序

这也是我需要公开道歉的地方。作为计算机科学家，当我们在一个“他们做错了”的领域工作时，我们偶尔会感到厌烦。地图制作就是一个例子。当我们到达现场时，地图总是记录在纬度上的位置，然后是经度;也就是说，南北轴线第一，东西通道第二。当你在屏幕上快速地画出它时，它看起来就像世界是横向的，因为在我看来坐标是“y/x”，你需要在画之前交换它们。

image.png

我们习惯于按照x/y的顺序工作，以至于我们最终会认为它应该是这样的——从那以后我们就一直在与地图制作者斗争。

所以如果你在“EPSG:4326”中看到一些数据，你不知道它是x/y顺序还是y/x顺序。

我们终于找到了一个替代方案;我们应该使用urn:ogc:def:crs:EPSG:6.6:4326而不是EPSG:4326。如果你曾经看到，你可以确定a)某人确实知道他们在做什么，b)数据被准确地记录在EPSG数据库定义的顺序。
变通方法
您可以在个案的基础上执行变通方法:

CRSAuthorityFactory   factory = CRS.getAuthorityFactory(true);
CoordinateReferenceSystem crs = factory.createCoordinateReferenceSystem("EPSG:4326");

或者你可以设置一个全局提示来强制GeoTools使用x/y顺序:

static void main(String args []){
   System.setProperty("org.geotools.referencing.forceXY", "true");

   ....
}

geotools学习（三）坐标系

坐标系

验证几何

导出重投影的shapefile

运行程序

导出重新投影的数据:

扩展尝试

Geometry

Point

Line

Polygon

为什么不使用Java Shape呢

坐标参考系

坐标系为我们定义了两个概念:

EPSG代码

轴顺序

相关文章

网友评论

延伸阅读

深度阅读

栏目导航

热点阅读

开源

geotools

geotools学习（三） 坐标系

坐标系

验证几何

导出重投影的shapefile

运行程序

导出重新投影的数据:

扩展尝试

Geometry

Point

Line

Polygon

为什么不使用Java Shape呢

坐标参考系

坐标系为我们定义了两个概念:

EPSG代码

轴顺序

相关文章

网友评论

延伸阅读

深度阅读

栏目导航

热点阅读

开源

geotools

geotools学习（三）坐标系