学前小故事

学前必备基础知识

WTK格式的Geomotry

GeoJSON格式的Geomotry

JTS(Java Topology Suite) Java拓扑套件

JTS(Java Topology Suite) 可视化界面

泰森多边形

向空间数据库插入数据

栅格

GDAL

GeoTools

QGIS

学前小故事

项目需求是跟用户当前位置判断是否在给定的地理位置范围内，符合位置限制才可以打卡，其中的位置范围是一个或多个不规则的多边形。如下图，判断用户是在清华还是北大。

图形获取区域坐标

因为项目前端使用微信小程序的wx.getLocation获取地理位置，为了坐标的一致性，后台选取区域范围采用了腾讯地图的地理位置服务，在应用工具->绘制几何图形里，提供了点、线、多边形和圆形可以方便的选取看这里

在官方提供的示例上稍加改动即可获取选定的位置坐标。

存储位置

取到坐标位置后，接着就是怎么存储？

开放地理空间联盟(OGC)是一个由 250多家公司，机构和大学组成的国际联盟，参与开发公开可用的空间解决方案，这些解决方案可用于管理空间数据的各种应用程序。OGC发布了地理信息的 OpenGIS®Implementation 标准，该规范可从 OGC 网站http://www.opengeospatial.org/standards/sfs获得。

为了遵循 OGC 规范，MySQL 将空间 extensions 实现为具有 Geometry Types 环境的 SQL 的子集，提供生成、存储、分析空间的功能。总之，MySQL可以满足我们的需求。

MySQL提供单个的存储类型 POINT、LINESTRING、POLYGON 对应几何图形点、线、多边形，GEOMETRY 可以存储三种中的任何一种。
同时拥有存储多种类型的能力， MULTIPOINT、MULTILINESTRING、MULTIPOLYGON、GEOMETRYCOLLECTION依次对应单个图形的复数。

回到项目中，我们用到的是 POLYGON ，建表语句 如下：

CREATE TABLE `polygon` (
  `id` int(10) unsigned NOT NULL AUTO_INCREMENT,
  `name` varchar(255) DEFAULT NULL,
  `polygon` polygon NOT NULL,
  PRIMARY KEY (`id`),
  SPATIAL KEY `d` (`polygon`)
) DEFAULT CHARSET=utf8;

插入数据

MySQL 支持将Well-Known 文本(WKT)格式和Well-Known 二进制(WKB)格式两种格式转换为object类型存储起来，我们使用更易于理解的WKT格式。

插入语句如下：

//GeomFromTextt()函数接受一个字符串，并且返回一个几何对象。
INSERT INTO `polygon` VALUES ('1', '清华大学', GeomFromText('POLYGON((
40.01169924229143 116.31565081888039,39.99304082299905 116.31616541796757,39.99343506780591 116.33297565023167,40.00237067000859 116.33743550702275,40.01340715321479 116.33057418815224,40.01169924229143 116.31565081888039))'));

INSERT INTO `polygon` VALUES ('2', '北京大学', GeomFromText('POLYGON((39.99711457525893 116.30450117461078,39.98673259872773 116.30535884106575,39.98673259872773 116.31702308311287,39.99963848242885 116.31598375134854,39.99711457525893 116.30450117461078))'));

需要注意的是腾讯地图返回的多边形的点不是闭合的，而polygon函数需要为了确定多边形是否闭合要求第一个点和最后一个点是一样的。如果不是闭合的polygon返回的结果将是NULL，插入语句就会执行失败。

如果几何满足诸如此(非穷举)列表中的条件，则它在语法上是 well-formed：

• 线串至少有两个点
• 多边形至少有一个环
• 多边形环关闭(第一个和最后一个点相同)
• 多边形环至少有 4 个点(最小多边形是一个三角形，第一个和最后一个点相同)
• 集合不为空(除了GeometryCollection)

查询判断

SELECT * FROM polygon WHERE
	MBRWithin (ST_GeomFromText('POINT(39.991333490218544 116.30964748487895)'), polygon);
# 在北京大学

SELECT * FROM polygon WHERE
	MBRWithin (ST_GeomFromText('POINT(39.988967560246685 116.3286905102832)'), polygon);
# 不在北大

细心的同学可能发现了这里的查询语句里用的是函数，在以往的SQL里如果存在查询字段上使用函数必然导致索引失效、全表扫描，但是在空间数据上不会，先看 EXPLAIN 语句和结果：

可见MySQL空间类型的数据同样可以建立索引，使用的关键词是 SPATIAL ，用法如下：

CREATE TABLE geom (g GEOMETRY NOT NULL);

CREATE SPATIAL INDEX g ON geom (g);

常用的空间计算函数

1、判断两点之间的距离
  ST_Distance(g1，g2)，返回g1和g2之间的距离。如果任一参数是NULL或空几何，则 return value 为NULL。
2、图形1是否完全包含图形2
  ST_Contains(g1，g2)，返回 1 或 0 以指示g1是否完全包含g2。还可以用ST_Within(g2,g1)达到相同的效果。
3、不相交
  ST_Disjoint(g1，g2)，返回 1 或 0 以指示g1是否在空间上与(不相交)g2不相交。

学前必备基础知识

1. 地理信息系统GIS

   地理信息系统（Geographic Information System或 Geo－Information system，GIS）有时又称为“地学信息系统”。它是一种特定的十分重要的空间信息系统。

   它是在计算机硬、软件系统支持下，对整个或部分地球表层（包括大气层）空间中的有关地理分布数据进行采集、储存、管理、运算、分析、显示和描述的技术系统。

2. ArcGIS平台

   ArcGIS产品线为用户提供一个可伸缩的，全面的GIS平台。ArcObjects包含了许多的可编程组件，从细粒度的对象（例如单个的几何对象）到粗粒度的对象（例如与现有ArcMap文档交互的地图对象）涉及面极广，这些对象为开发者集成了全面的GIS功能。

3. 矢量化数据

   矢量数据是在直角坐标中，用x、y坐标表示地图图形或地理实体的位置和形状的数据。矢量数据一般通过记录坐标的方式来尽可能地将地理实体的空间位置表现得准确无误。在arcgis中数据矢量化就是将栅格数据转为矢量数据。

   arcgis作为地理信息行业的功能之一数据矢量化的详细步骤：

   • 添加数据按钮打开栅格数据和待矢量用的GDB数据或者shp数据或者mdb数据（用于存矢量化后的独迅称数据）
   • 编辑矢量数据，右键点击数据——编辑要素——开始编辑
   • 使用鼠标根据栅格影像中的地物形状进行描绘，一般就是在拐点处添加节点（点击鼠标左键即添加节点），描绘完成后，双击即可
   • 再对刚刚完成的数据进行赋属性，也就是在属性表中输入描绘的类型，如这里是“房子”。
   • 再使用鼠标对下一个地物进行描绘，如图中的公路。以此循环对影像中各种地物进行描绘，这便是矢量化的过程

4. Geometry数据类型

   Geometry是一种空间几何数据类型，常用于描述空间几何信息，例如坐标点、线、面、三维信息等。

   也就是说，GIS一般使用Geometry数据类型来存储及展示地理信息。

   常见的支持Geometry的数据库有Oracle、SqlServer、Mysql、PostgreSQL

5. 实际运用特殊性

   数据库中有geometry类型，但是很多语言中没有直接的geometry类型。
   以java为例，如果想要存储geometry字段，需要增加SET @g1 = geomFromText(geometry) ，也就是将字符串转换成几何对象存储；

   同样，从数据库中直接取出的geometry类型字段，java无法使用，需要在数据库中查询的时候，做个几何类型转字符串类型的处理st_AsText(geom) as geom

6. Geometry之间的关系

   关系	说明
   相等（Equals）	几何形状拓扑上相等
   脱节（Disjoint）	几何形状没有共有的点
   相交（Intersects）	几何形状至少有一个共有点（区别于脱节）
   接触（Touches）	几何形状有至少一个公共的边界点，但是没有内部点
   交叉（Crosses）	几何形状共享一些但不是所有的内部点
   内含（Within）	几何形状A的线都在几何形状B内部
   包含（Contains）	几何形状B的线都在几何形状A的内部（区别于内含）
   重叠（Overlaps）	几何形状共享一部分但不是所有的公共点，而且相交处有他们自己相同的区域

7. 关系分析

   分析	说明
   缓冲区分析（Buffer）	包含所有的点在一个指定距离内的多边形和多多边形
   凸壳分析（ConvexHull）	包含几何形体的所有点的最小凸壳多边形（外包多边形）
   交叉分析（Intersection）	A∩B 交叉操作就是多边形AB中所有共同点的集合
   联合分析（Union）	AUB AB的联合操作就是AB所有点的集合
   差异分析（Difference）	(A-A∩B) AB形状的差异分析就是A里有B里没有的所有点的集合
   对称差异分析（SymDifference）	(AUB-A∩B) AB形状的对称差异分析就是位于A中或者B中但不同时在AB中的所有点的集合

新增空间数据

当我们有了数据库，也有了arcgis提供的服务，现在该如何新增空间数据呢？

先通过WKT了解下空间数据Geometry大概长什么样

WKT，是一种文本标记语言，用于表示矢量几何对象、空间参照系统及空间参照系统之间的转换。
它的二进制表示方式，亦即WKB(well-known binary)则胜于在传输和在数据库中存储相同的信息。该格式由开放地理空间联盟(OGC)制定。

WKT可以表示的几何对象包括：点，线，多边形，TIN（不规则三角网）及多面体。可以通过几何集合的方式来表示不同维度的几何对象。几何物体的坐标可以是2D(x,y),3D(x,y,z),4D(x,y,z,m),加上一个属于线性参照系统的m值。

以下为几何WKT字串样例：
◈ POINT(6 10)
◈ LINESTRING(3 4,10 50,20 25)
◈ POLYGON((1 1,5 1,5 5,1 5,1 1),(2 2,2 3,3 3,3 2,2 2))
◈ MULTIPOINT(3.5 5.6, 4.8 10.5)
◈ MULTILINESTRING((3 4,10 50,20 25),(-5 -8,-10 -8,-15 -4))
◈ MULTIPOLYGON(((1 1,5 1,5 5,1 5,1 1),(2 2,2 3,3 3,3 2,2 2)),((6 3,9 2,9 4,6 3)))
◈ GEOMETRYCOLLECTION(POINT(4 6),LINESTRING(4 6,7 10))
◈ POINT ZM (1 1 5 60)
◈ POINT M (1 1 80)
◈ POINT EMPTY
◈ MULTIPOLYGON EMPTY

到这里我们清楚了，所谓空间数据就是一个或一组坐标，这个坐标或坐标组有类型（POINT点LINESTRING线POLYGON面），通过这些坐标，GIS系统可以完整地定位查询绘制这些坐标信息

Java可以将WKT格式的Geomotry转换成GeoJSON

1. Meven添加依赖

   <!-- 引入json处理包 -->
   <dependency>
       <groupId>com.alibaba</groupId>
       <artifactId>fastjson</artifactId>
       <version>1.2.47</version>
   </dependency>
   
   <!-- jts处理Geometry -->
   <dependency>
       <groupId>com.vividsolutions</groupId>
       <artifactId>jts</artifactId>
       <version>1.13</version>
   </dependency>
   
   <!-- https://mvnrepository.com/artifact/org.geotools/gt-geojson -->
   <dependency>
       <groupId>org.geotools</groupId>
       <artifactId>gt-geojson</artifactId>
       <version>9.3</version>
   </dependency>
   
   <dependency>
       <groupId>com.googlecode.json-simple</groupId>
       <artifactId>json-simple</artifactId>
       <version>1.1.1</version>
   </dependency>
   

   如添加依赖遇到not fount，可手动配置maven repository解决

   <repositories>
       <repository>
           <id>OSGeo Repository</id>
           <url>http://download.osgeo.org/webdav/geotools/</url>
       </repository>
   </repositories>
   

2. 解析方法源码

   package com.bret.utils;
   
   import com.vividsolutions.jts.geom.Geometry;
   import com.vividsolutions.jts.io.WKTReader;
   import org.geotools.geojson.geom.GeometryJSON;
   
   import java.io.StringWriter;
   
   public class WKTUtil {
   
       /**
        * 由wkt格式的geometry生成geojson
        * @param wkt
        * @return
        */
       public static String wktToJson(String wkt) {
           String json = null;
           try {
               WKTReader reader = new WKTReader();
               Geometry geometry = reader.read(wkt);
               StringWriter writer = new StringWriter();
               GeometryJSON g = new GeometryJSON(20);
               g.write(geometry, writer);
               json = writer.toString();
           } catch (Exception e) {
               e.printStackTrace();
           }
           return json;
       }
   
   }
   

   POINT(6 10) --> {“type”:“Point”,“coordinates”:[6,10]}
   LINESTRING(3 4,10 50,20 25) --> {“type”:“LineString”,“coordinates”:[[3,4],[10,50],[20,25]]}

GeoJSON格式的Geomotry

GeoJSON是一种对各种地理数据结构进行编码的格式。GeoJSON对象可以表示几何、特征或者特征集合。
GeoJSON支持下面几何类型：点、线、面、多点、多线、多面和几何集合。GeoJSON里的特征包含一个几何对象和其他属性，特征集合表示一系列特征。

GeoJSON 可以转 shp格式，GeoTools提供封装方法，详情可自行查询。(GeoJSON --> .shp --> .nc)

一个完整的GeoJSON数据结构总是一个（JSON术语里的）对象。在GeoJSON里，对象由名/值对–也称作成员的集合组成。对每个成员来说，名字总是字符串。成员的值要么是字符串、数字、对象、数组，要么是下面文本常量中的一个：“true”,“false"和"null”。数组是由值是上面所说的元素组成。

GeoJSON特征集合举例：

{
  //type是由GeoJSON对象的类型所确定的字符串，类型有："Point", "MultiPoint", "LineString", "MultiLineString", "Polygon", "MultiPolygon",   	     "GeometryCollection", "Feature", "FeatureCollection"
  "type": "FeatureCollection", 
  //特征对象的结合数组
  "features": [
    { "type": "Feature",
      "geometry": {"type": "Point", "coordinates": [102.0, 0.5]},
      "properties": {"prop0": "value0"}
      },
    { "type": "Feature",
      "geometry": {
        "type": "LineString",
        "coordinates": [
          [102.0, 0.0], [103.0, 1.0], [104.0, 0.0], [105.0, 1.0]
          ]
        },
      "properties": {
        "prop0": "value0",
        "prop1": 0.0
        }
      },
    { "type": "Feature",
       "geometry": {
         "type": "Polygon",
         "coordinates": [
           [ [100.0, 0.0], [101.0, 0.0], [101.0, 1.0],
             [100.0, 1.0], [100.0, 0.0] ]
           ]
       },
       "properties": {
         "prop0": "value0",
         "prop1": {"this": "that"}
         }
       }
     ]
   }

• featureCollection（特征集合）：表示包含feature（空间有界实体）对象
• features表示所包含一个个feature对象
• type：表示当前层的类型
• properties：表示当前层的信息，比如name = 上城区
• geometry：几何的值

JTS(Java Topology Suite) Java拓扑套件

wkt文本仅仅是一个字符串而已，直接将坐标点拼接成符合WKT格式的字符串不就可以了吗？道理是这个道理，要做好可就难了。难点如下:

• 拼接工作量巨大
• 拼接过程容易出错
• 拼接的结果不一定合法可用

对此我们需要一套JAVA API对数据进行处理，该API能够方便的创建Geometry对象，进行地理信息的绘制、创建、验证等等功能，市面上常见的GeometryApi有

① Esri/geometry-api-java
② locationtech/jts（推荐）

Esri是Arcgis官方提供的javaSDK，可惜功能不多，甚至不能提供基本的空间计算功能。jts功能较为齐全，资料也相对丰富一点。这里我们就已JTS为重点进行讲解。

JTS(Java Topology Suite) Java拓扑套件，是Java的处理地理数据的API。它的API文档地址：https://locationtech.github.io/jts/javadoc/

JTS包结构如下：

• 系（linearref包）
• 计算交点（noding包）
• 几何图形操作（operation包）
• 平面图（planargraph包）
• 多边形化（polygnize包）
• 精度（precision）
• 工具（util包）

JTS提供的空间数据类型如下：

数据类型	内容
Point	点
MultiPoint	多点
LineString	线
LinearRing	封闭的线条
MultiLineString	多条线条
Polygon	多边形
MultiPolygon	多个多边形
GeometryCollection	包括点，线，面

JTS提供的空间拓扑关系如下：

空间关系	中文名称	OGC标准	解释
Contains	包含	是	一个几何图形的内容完全包含了另一个几个图形的内部和边界
CoveredBy	覆盖	否	一个几何图形被另一个几何图形所包含，并且它们的边界相交。Point和MultiPoint不支持此空间关系,因为它们没有边界
Crosses	交叉	是	一个几何图形的内部和另一个几何图形的边界和内部相交，但它们的边界不相交.（两条线相交于一点，一条线和一个面相交）
Disjoint	分离	是	两个几何图形的边界和内部不相交
EnvelopeIntersects	封套相交	否	两个几个图形的外接矩形相交
Equal	相等	是	两个几何图形具有相同的边界和内部
Inside	内部	否	一个几何图形在另一个几何图形的内部，但是和它的边界不接触
Intersects	相交	是	两个几何图形没有分离（Non-DisJoint）
Overlaps	重叠	是	两个几何图形的边界和内部相交（Intersect）
Touch	接触	是	两个几何图形的边界相交,但是内部不相交
Within	包含于	是	一个几何图形的内部和边界完全在另一个几何图形的内部

JTS提供的几何叠加操作如下：

空间关系	中文名称	OGC标准	解释
Buffer	缓冲区分析	是	几何图形指定距离（扩大或者缩小）后的几何图形
ConvexHull	凸壳分析	是	包含几何形体的所有点的最小凸壳多边形（外包多边形）
Intersection	交叉分析	是	交叉操作就是多边形AB中所有共同点的集合
Union	联合分析	是	AB的联合操作就是AB所有点的集合
Difference	差异分析	是	AB形状的差异分析就是A里有B里没有的所有点的集合
SymDifference	对称差异分析	是	AB形状的对称差异分析就是位于A中或者B中但不同时在AB中的所有点的集合

我们在使用JTS前首先需要引入JTS的Maven坐标：

<!-- https://mvnrepository.com/artifact/org.locationtech.jts/jts-core -->
<dependency>
    <groupId>org.locationtech.jts</groupId>
    <artifactId>jts-core</artifactId>
    <version>1.18.2</version>
</dependency>

使用JTS前先了解下Geometry类常用工具方法

1. GeometryFactory

   GeometryFactory是geometry的工厂类，提供了Geometry子类的创建方法的统一入口。在创建Geometry的时候建议用GeometryFactory去创建，不要单独自己去new一个子类。

   数据结构:
   private PrecisionModel precisionModel : 精度模式，默认为FLOAT
   private CoordinateSequenceFactory coordinateSequenceFactory : 坐标序列工厂类，创CoordinateList

2. PrecisionModel

   创建Geometry的时候通过PrecisionModel可以控制geometry坐标的精度。
   三种模式:
   FLOATING: 等同于Java的双精度double类型
   FLOATING_SINGLE: 等同于Java的单精度float类型
   FIXED: 固定精度模式(固定小数点后几位)。通过public PrecisionModel(double scale)构造方法构造。
   FIXED去固定精度并不会改变geometry里的点的真实坐标，只会去改变toString出来的wkt格式的坐标精度。慎用。

3. LineString

   数据结构：CoordinateSequence coordinates，即一堆点的集合。

   method	description
   boolean isClosed()	判断LineString是否闭合,第一个点和最后一个点是否相等
   boolean isRing()	在判断isClosed()的同时判断是否是简单图形
   LineString copyInternal()	深拷贝一个LineString

4. LinearRing

5. Polygon

   数据结构：
   LinearRing shell : 表示外轮廓。一个Polygon只能有一个shell。
   LinearRing[] holes: 表示洞。一个Polygon可有多个hole。

   method	description
   Coordinate[] getCoordinates()	获取Polygon所有的点，polygon中洞的点会排在外轮廓的点之后
   boolean isRectangle()	判断polygon是不是矩形
   LinearRing getExteriorRing()	获取polygon外轮廓
   LinearRing getInteriorRingN(int n)	获取某一个洞
   void apply(CoordinateFilter filter)	直接修改一个polygon的值(比如精度修改，放大缩小之类的操作)
   Polygon copyInternal()	等同于深拷贝
   Geometry convexHull()	拿到Polygon的凸包
   void normalize()	将polygon格式化为逆时针构建,一般equal之前会先normalize一下
   Geometry normalize()	同上，只是返回结果为一个Polygon

6. GeometryCollection

   数据结构：
   Geometry[] geometries: geometry的集合。子类有MultiPoint MultiLineString MultiPolygon.
   实现了所有的Geometry的抽象方法，只是改成了对集合整体的操作

现在展示一下点、线、多边形的JTS创建过程：

1. 创建点

   
   private GeometryFactory geometryFactory = JTSFactoryFinder.getGeometryFactory( null );
   
   /**
    * 创建点
    * @return
    */
   public Point createPoint(){
   
   	Coordinate coord = new Coordinate(109.013388, 32.715519);
   	Point point = new GeometryFactory().createPoint(coord);
   	return point;
   }
   
   /**
    * 创建点（通过WKT）
    * @return
    */
   public Point createPointByWKT() throws ParseException{
   
   	WKTReader reader = new WKTReader( geometryFactory );
   	Point point = (Point) reader.read("POINT (109.013388 32.715519)");
   	return point;
   }
   

2. 创建线

   
   private GeometryFactory geometryFactory = JTSFactoryFinder.getGeometryFactory( null );
   
   /**
    * 创建线
    * @return
    */
   public LineString createLine(){
   
   	Coordinate[] coords  = new Coordinate[] {new Coordinate(2, 2), new Coordinate(2, 2)};
   	LineString line = geometryFactory.createLineString(coords);
   	return line;
   }
   
   /**
    * 创建线（通过WKT）
    * @return
    */
   public LineString createLineByWKT() throws ParseException{
   
   	WKTReader reader = new WKTReader( geometryFactory );
   	LineString line = (LineString) reader.read("LINESTRING(0 0, 2 0)");
   	return line;
   }	
   

3. 创建多条线

   
   private GeometryFactory geometryFactory = JTSFactoryFinder.getGeometryFactory( null );
   
   /**
    * 创建多条线
    * @return
    */
   public MultiLineString createMLine(){
   
       Coordinate[] coords1  = new Coordinate[] {new Coordinate(2, 2), new Coordinate(2, 2)};
       LineString line1 = geometryFactory.createLineString(coords1);
       
       Coordinate[] coords2  = new Coordinate[] {new Coordinate(2, 2), new Coordinate(2, 2)};
       LineString line2 = geometryFactory.createLineString(coords2);
       
       LineString[] lineStrings = new LineString[2];
       lineStrings[0]= line1;
       lineStrings[1] = line2;
       
       MultiLineString ms = geometryFactory.createMultiLineString(lineStrings);
       return ms;
   }
   
   
   /**
    * 创建多条线（通过WKT）
    * @return
    * @throws ParseException
    */
   public MultiLineString createMLineByWKT()throws ParseException{
   
       WKTReader reader = new WKTReader( geometryFactory );
       MultiLineString line = (MultiLineString) reader.read("MULTILINESTRING((0 0, 2 0),(1 1,2 2))");
       return line;
   }
   

4. 创建多边形

   
   private GeometryFactory geometryFactory = JTSFactoryFinder.getGeometryFactory( null );
   
   /**
    * 创建多边形（通过WKT）
    * @return
    * @throws ParseException
    */
   public Polygon createPolygonByWKT() throws ParseException{
       WKTReader reader = new WKTReader( geometryFactory );
       Polygon polygon = (Polygon) reader.read("POLYGON((20 10, 30 0, 40 10, 30 20, 20 10))");
       return polygon;
   }
   
   /**
    * 创建多个多边形（通过WKT）
    * @return
    * @throws ParseException
    */
   public MultiPolygon createMulPolygonByWKT() throws ParseException{
       WKTReader reader = new WKTReader( geometryFactory );
       MultiPolygon mpolygon = (MultiPolygon) reader.read("MULTIPOLYGON(((40 10, 30 0, 40 10, 30 20, 40 10),(30 10, 30 0, 40 10, 30 20, 30 10)))");
       return mpolygon;
   }
   
   

现在展示一下几何图形之间的拓扑关系：

1. 相等（Eequals）

   /**
    *  两个几何对象是否是重叠的
    * @return
    * @throws ParseException
    */
   public boolean equalsGeo() throws ParseException{
       WKTReader reader = new WKTReader( geometryFactory );
       LineString geometry1 = (LineString) reader.read("LINESTRING(0 0, 2 0, 5 0)");
       LineString geometry2 = (LineString) reader.read("LINESTRING(5 0, 0 0)");
       return geometry1.equals(geometry2);//true
   }
   

2. 脱节（Disjoint）

   /**
    * 几何对象没有交点(相邻)
    * @return
    * @throws ParseException
    */
   public boolean disjointGeo() throws ParseException{
       WKTReader reader = new WKTReader( geometryFactory );
       LineString geometry1 = (LineString) reader.read("LINESTRING(0 0, 2 0, 5 0)");
       LineString geometry2 = (LineString) reader.read("LINESTRING(0 1, 0 2)");
       return geometry1.disjoint(geometry2);
   }
   
   

3. 相交（Intersects）

   /**
    * 至少一个公共点(相交)
    * @return
    * @throws ParseException
    */
   public boolean intersectsGeo() throws ParseException{
       WKTReader reader = new WKTReader( geometryFactory );
       LineString geometry1 = (LineString) reader.read("LINESTRING(0 0, 2 0, 5 0)");
       LineString geometry2 = (LineString) reader.read("LINESTRING(0 0, 0 2)");
       Geometry interPoint = geometry1.intersection(geometry2);//相交点
       System.out.println(interPoint.toText());//输出 POINT (0 0)
       return geometry1.intersects(geometry2);
   }
   

4. 内含（within）

   /**
    * 判断以x,y为坐标的点point(x,y)是否在geometry表示的Polygon中
    * @param x
    * @param y
    * @param geometry wkt格式
    * @return
    */
   public boolean withinGeo(double x,double y,String geometry) throws ParseException {
   
       Coordinate coord = new Coordinate(x,y);
       Point point = geometryFactory.createPoint( coord );
   
       WKTReader reader = new WKTReader( geometryFactory );
       Polygon polygon = (Polygon) reader.read(geometry);
       return point.within(polygon);
   }
   

5. 包含（Contains）

   public boolean containsGeo() throws ParseException {
   
   	Point point = new GeometryFactory().createPoint(new Coordinate(1, 1));
   	return geometry1.equals(geometry2);
   
   }
   

现在展示一下几何图形的叠加操作：

1. 缓冲区分析（Buffer）

   一条linestring经过buffer之后某种意义上变成一个ploygon，下面蓝色的线经过buffer之后成为了这个样子，如此就可以判断两条在一定阈值内重合的线，实质上是扩大了一条线的范围，然后再求buffer这个ploygon与另外一条linestring有无交集。

   

   /**
    * 返回a指定距离内的多边形和多多边形
    *
    * @param a
    * @param distance
    * @return
    */
   public Geometry bufferGeo(Geometry a, double distance) {
       return a.buffer(distance);
   }
   

2. 交叉分析（intersection）

   计算出预测线与实际线再阈值为1的情况下重合的部分：
   

   /**
    * 创建g1-预测值
    *
    * @param a
    * @param b
    * @return
    */
   public static void main(String[] args) throws ParseException {
   
   	Geometry g1 = new WKTReader().read("LINESTRING (0 0, 2 3, 3 2.5, 5 4, 6 3)");
   	System.out.println("Geometry 1: " + g1);
   	
   	// 创建g2-实际值
   	Geometry g2 = new WKTReader().read("LINESTRING (1 1, 2 2, 3 1, 5 3, 6 2, 8 3)");
   	System.out.println("Geometry 2: " + g2);
   	
   	Geometry buffer = g1.buffer(1);
   	// 重合的部分
   	Geometry g1_g2 = buffer.intersection(g2);
   	System.out.println("G1 intersection G2: " + g1_g2);
   }
   

   

3. 联合分析（union）

   /**
    * 几何对象合并
    *
    * @param a
    * @param b
    * @return
    */
   public Geometry unionGeo(Geometry a, Geometry b) {
       return a.union(b);
   }
   

4. 差异分析（difference）

   根据示例2中内容拓展

   /**
    * 创建g1-预测值
    *
    * @param a
    * @param b
    * @return
    */
   public Geometry differenceGeo(Geometry a, Geometry b) {
   
   	Geometry g1 = new WKTReader().read("LINESTRING (0 0, 2 3, 3 2.5, 5 4, 6 3)");
   	System.out.println("Geometry 1: " + g1);
   	
   	// 创建g2-实际值
   	Geometry g2 = new WKTReader().read("LINESTRING (1 1, 2 2, 3 1, 5 3, 6 2, 8 3)");
   	System.out.println("Geometry 2: " + g2);
   	
   	Geometry buffer = g1.buffer(1);
   	// 重合的部分
   	Geometry g1_g2 = buffer.intersection(g2);
   	System.out.println("G1 intersection G2: " + g1_g2);
   	
   	// g2- g1_g2
   	Geometry difference_addLine = g2.difference(g1_g2);
   	System.out.println("G2 difference g1_g2: " + difference_addLine);
   }
   

   

5. 支持接口

   • Coordinate

     Coordinate（坐标）是用来存储坐标的轻便的类。它不同于点，点是Geometry的子类。不像模范Point的对象（包含额外的信息，例如一个信包，一个精确度模型和空间参考系统信息），Coordinate只包含纵坐标值和存取方法。

   • Envelope(矩形)

     一个具体的类，包含一个最大和最小的x 值和y 值。

   • GeometryFactory

     GeometryFactory提供一系列的有效方法用来构造来自Coordinate类的Geometry对象。支持接口

JTS(Java Topology Suite) 可视化界面

JTS提供了一个界面工具JTS TestBuilder，可以在上面绘制图形，验证各种算法。

链接：https://pan.baidu.com/s/1cui-yyGeHRG6bW6EovVb9g
提取码：pmr7

下载后将压缩包解压，在jts-1.14文件夹下的bin文件夹中，双击testbuilder.bat，既可启动JTS TestBuilder。

侧边工具栏里的内容概览：

• Valid/Mark用于判断几何图形是否合法简单
• Predicates获取DE-9IM模型结果；
• Scalar Function查验几何图形的空间关系；
• Geometry Functions对几何对象进行仿射、缓冲等操作，并获取结果

底部工具栏里的内容概览：

• Cases：测试用例界面
• Input：用来输入图形参数的，目前只有A和B两个参数，默认测试函数都是用参数A
• Result：计算的结果
• Inspect：可以参看geometry的具体vetex数据，方向(在Input的Inspect按钮点击进行查看当前高亮的编辑框的对象)
• Layers：可以控制输入图形和结果图层的显示，方便查看结果

JTS TestBuilder的使用方法

1. 判断两个面是否相交

   点击工具条上的Draw Polygon

   

   进入侧边栏，点击input，点击A输入框，在图板上绘制图形，双击结束绘制

   
   点击input的B输入框，切换到绘制图形B，绘制完成后，双击结束

   
   在Scalar Functions中，选择Intersects，点击compute，下面控制台中输出结果

   

泰森多边形

泰森多边形的构造

Geotools生成泰森多边形

1. 泰森多边形的构造

   假如我在一片空地上，空地上有两个厕所

   

   此时我突然内急，为了赶时间我肯定要去离自己近的那个，我们把厕所视作控制点，然后用个控制点连线的中垂线把整个区域一分为二

   
   然后根据自己所处的位置选择离自己更近的那一个，当再增加一个控制点，上面的问题会变得稍微复杂一些

   
   我们通过三条中垂线把区域分成3份，每块区域内都有对应的最佳选择，当控制点发生变化，3块区域也会随之调整变化。

   这样的图形被称为3控制点沃罗诺伊图，图中的每一个单元称为沃罗诺伊单元，或者也可以被称为泰森多边形

   沃罗诺伊图满足下面三个特征：
   • 每个泰森多边形内仅含有一个控制点
   • 泰森多边形内的点到相应控制点的距离最近
   • 位于泰森多边形边上的点到其两边的控制点的距离相等

   沃罗诺伊图有非常多的实际应用案例：

   • 基于全球机场的沃罗诺伊图，每座机场都是一个控制点，把地球分割成多个泰森多边形，这样飞机能够时刻知晓离自己最近的机场

   • 基站的建设

2. Geotools生成泰森多边形

   首先导入相关依赖：

   <properties>
       <geotools.version>19.4</geotools.version>
   </properties>
   
   <dependencies>
   	<!-- Geotools -->
   	<dependency>
     		<groupId>org.geotools</groupId>
   		<artifactId>gt-shapefile</artifactId>
   		<version>${geotools.version}</version>
       </dependency>
       <dependency>
   		<groupId>org.geotools</groupId>
   		<artifactId>gt-geojson</artifactId>
   		<version>${geotools.version}</version>
       </dependency>
       <dependency>
   		<groupId>org.geotools</groupId>
   		<artifactId>gt-main</artifactId>
   		<version>${geotools.version}</version>
   	</dependency>
       <dependency>
   		<groupId>org.geotools</groupId>
   		<artifactId>gt-epsg-hsql</artifactId>
   		<version>${geotools.version}</version>
   	</dependency>
   	<dependency>
   	    <groupId>org.geotools</groupId>
   	    <artifactId>gt-geotiff</artifactId>
   	    <version>${geotools.version}</version>
   	</dependency>
   	<dependency>
   	    <groupId>org.geotools</groupId>
   	    <artifactId>gt-process-raster</artifactId>
   	    <version>${geotools.version}</version>
   	</dependency>
   	<!-- https://mvnrepository.com/artifact/org.geotools/gt-image -->
   	<dependency>
   	    <groupId>org.geotools</groupId>
   	    <artifactId>gt-image</artifactId>
   	    <version>${geotools.version}</version>
   	</dependency>
   </dependencies>
   

   创建测试点：

   

   int xmin = 0, xmax=180;
   int ymin = 0, ymax=90;
   Random random = new Random();
   List<Geometry> geomsPoints = new ArrayList<Geometry>();
   for(int i=0;i<100;i++){
   	//生成0-100之间随机数（制作模拟点）
   	int x = random.nextInt(xmax)%(xmax-xmin+1) + xmin,
           y = random.nextInt(ymax)%(ymax-ymin+1) + ymin;
       Coordinate coord = new Coordinate(x,y,i);
       coords.add(coord);
       clipEnvelpoe.expandToInclude(coord);
       geomsPoints.add(new GeometryFactory().createPoint(coord));
   }
   

   生成泰森多边形：

   

   voronoiDiagramBuilder.setSites(coords);
   voronoiDiagramBuilder.setClipEnvelope(clipEnvelpoe);
   Geometry geom = voronoiDiagramBuilder.getDiagram(JTSFactoryFinder.getGeometryFactory());
   List<Geometry> geoms = new ArrayList<Geometry>();
   for(int i=0;i<geom.getNumGeometries();i++){
       geoms.add(geom.getGeometryN(i));
   }
   

   完整代码如下：

   package com.lzugis.geotools;
   
   import java.io.File;
   import java.io.Serializable;
   import java.nio.charset.Charset;
   import java.util.ArrayList;
   import java.util.HashMap;
   import java.util.List;
   import java.util.Map;
   import java.util.Random;
   
   import org.geotools.data.FeatureWriter;
   import org.geotools.data.Transaction;
   import org.geotools.data.shapefile.ShapefileDataStore;
   import org.geotools.data.shapefile.ShapefileDataStoreFactory;
   import org.geotools.feature.simple.SimpleFeatureTypeBuilder;
   import org.geotools.geometry.jts.JTSFactoryFinder;
   import org.geotools.referencing.crs.DefaultGeographicCRS;
   import org.opengis.feature.simple.SimpleFeature;
   import org.opengis.feature.simple.SimpleFeatureType;
   
   import com.vividsolutions.jts.geom.Coordinate;
   import com.vividsolutions.jts.geom.Envelope;
   import com.vividsolutions.jts.geom.Geometry;
   import com.vividsolutions.jts.geom.GeometryFactory;
   import com.vividsolutions.jts.geom.Point;
   import com.vividsolutions.jts.geom.Polygon;
   import com.vividsolutions.jts.triangulate.VoronoiDiagramBuilder;
   
   public class TsdbxTest {
   	static TsdbxTest tsdbx = new TsdbxTest();
   	public void voronoiTest(){
   	    VoronoiDiagramBuilder voronoiDiagramBuilder = new VoronoiDiagramBuilder();
   	    List<Coordinate> coords = new ArrayList<Coordinate>();
   	    Envelope clipEnvelpoe = new Envelope();
   	    int xmin = 0, xmax=180;
           int ymin = 0, ymax=90;
           Random random = new Random();
           List<Geometry> geomsPoints = new ArrayList<Geometry>();
   	    for(int i=0;i<100;i++){
   	    	int x = random.nextInt(xmax)%(xmax-xmin+1) + xmin,
   	            y = random.nextInt(ymax)%(ymax-ymin+1) + ymin;
   	        Coordinate coord = new Coordinate(x,y,i);
   	        coords.add(coord);
   	        clipEnvelpoe.expandToInclude(coord);
   	        geomsPoints.add(new GeometryFactory().createPoint(coord));
   	    }
   	    String pointpath = "d:/data/tsdbxpt.shp";
   	    tsdbx.writeShape(pointpath,"Point", geomsPoints);
   
   	    voronoiDiagramBuilder.setSites(coords);
   	    voronoiDiagramBuilder.setClipEnvelope(clipEnvelpoe);
   	    Geometry geom = voronoiDiagramBuilder.getDiagram(JTSFactoryFinder.getGeometryFactory());
   	    List<Geometry> geoms = new ArrayList<Geometry>();
   	    for(int i=0;i<geom.getNumGeometries();i++){
   	        geoms.add(geom.getGeometryN(i));
   	    }
   	    
   	    String polygonpath = "d:/data/tsdbx.shp";
   	    tsdbx.writeShape(polygonpath,"Polygon", geoms);
   	}
   	/**
   	 * 
   	 * @param filepath
   	 * @param geoType
   	 * @param geoms
   	 */
   	public void writeShape(String filepath, String geoType, List<Geometry> geoms) {
   		try {
   			//创建shape文件对象
   			File file = new File(filepath);
   			Map<String, Serializable> params = new HashMap<String, Serializable>();
   			params.put( ShapefileDataStoreFactory.URLP.key, file.toURI().toURL() );
   			ShapefileDataStore ds = (ShapefileDataStore) new ShapefileDataStoreFactory().createNewDataStore(params);
   			//定义图形信息和属性信息
   			SimpleFeatureTypeBuilder tb = new SimpleFeatureTypeBuilder();
   			tb.setCRS(DefaultGeographicCRS.WGS84);
   			tb.setName("shapefile");
   			if(geoType=="Point"){
   				tb.add("the_geom", Point.class);
   			}
   			else{
   				tb.add("the_geom", Polygon.class);
   			}
   			
   			ds.createSchema(tb.buildFeatureType());
   			//设置编码
               Charset charset = Charset.forName("GBK");
               ds.setCharset(charset);
   			//设置Writer
   			FeatureWriter<SimpleFeatureType, SimpleFeature> writer = ds.getFeatureWriter(ds.getTypeNames()[0], Transaction.AUTO_COMMIT);
   			for(int i=0,len=geoms.size();i<len;i++){
   				//写下一条
   				SimpleFeature feature = writer.next();
   				Geometry geom = geoms.get(i);
   				feature.setAttribute("the_geom", geom);
   			}
   			writer.write();
   			writer.close();
   			ds.dispose();
   		} 
   		catch (Exception e) {
   			e.printStackTrace();
   		}
   	}
   	
   	public static void main(String[] args){
   		long start = System.currentTimeMillis();
   		tsdbx.voronoiTest();
   		System.out.println("共耗时"+(System.currentTimeMillis() - start)+"ms");
   	}
   }
   

   在geotools中是内置了泰森多边形的算法的，主要依赖于 VoronoiDiagramBuilder 类 利用如下的方式就完成了泰森多边形的计算：

   VoronoiDiagramBuilder voronoiDiagramBuilder = new VoronoiDiagramBuilder();
   		List<Coordinate> coords = new ArrayList<Coordinate>();
           //这是读取了包含elevation字段的点数据
   		SimpleFeatureIterator simpleFeatureIterator = featureCollection.features();
   		while(simpleFeatureIterator.hasNext()) {
   		    SimpleFeature simpleFeauture = simpleFeatureIterator.next();
   		    Point point = (Point)simpleFeauture.getDefaultGeometry();
   		    double elevation = Double.parseDouble(simpleFeauture.getProperty("elevation").getValue().toString());
   			Coordinate coord = new Coordinate(point.getX(),point.getY(),elevation);
   			coords.add(coord);
   		}
           simpleFeatureIterator.close();
           //输入坐标点
           voronoiDiagramBuilder.setSites(coords);
           //设置泰森多边形的边界（一个com.vividsolutions.jts.geom.Polygon）
           voronoiDiagramBuilder.setClipEnvelope(boundary.getEnvelopeInternal());
           //输出成果
           Geometry geom = voronoiDiagramBuilder.getDiagram(gf);
   

向空间数据库插入数据

数据库如何新增这些空间数据呢？

先了解一下SQLSERVER空间数据的插入语句长啥样

--GEOM是类型为Geometry的字段--
--我们向该字段新增了一条3D的多边形数据--
--geometry :: STGeomFromText () 是由SQLSERVER提供的函数，它能将WKT文本转换为数据库geometry类型的数据--
INSERT INTO [dbo].[TEST_GEO_TABLE] ( [GEOM] )
VALUES
    ( geometry :: STGeomFromText ( 
    'POLYGON ((
        113.507259000000005 22.24814946 8, 
        113.507188600000006 22.248088559999999 9, 
        113.507117399999998 22.24802743 10, 
        113.507046099999997 22.24796624 11, 
        113.507017300000001 22.247888209999999 12
        ))',4326 )
    );

也就是说，将坐标转化为WKT文本，我们就可以插入空间数据。接下来我们要考虑的是如何产生WKT文本

栅格

ArcGIS栅格数据集通过将世界分割成在格网上布局的离散方形或矩形像元来表示地理要素。每个像元都具有一个值，用于表示该位置的某个特征，例如温度、高程或光谱值。可以抽象为数学上的“矩阵”，这一点与数字图像类似，与其不同的是栅格数据往往带有空间信息（地理坐标、投影坐标）。

像元即栅格所表示的内容的详细程度（要素/现象）通常取决于像元（像素）大小或空间分辨率。像元必须足够小，这样才可以捕获到所需的详细信息；而像元又必须足够大，这样才可以提高计算机存储和分析的执行效率。栅格可以使用更小的像元大小在要素的范围内表示更多的特征、更小的要素或更详细的内容。不过，更多通常未必更好。像元大小如果较小，则在表示整个表面时会造成栅格数据集较大；因此，会需要更大的储存空间，而且通常会使处理时间更长。

栅格数据集常用于表示和管理影像、数字高程模型及许多其他现象。通常，栅格是用于表示点、线和多边形要素的一种方法。在下面的示例中，可以看到如何将一系列多边形表示为栅格数据集。

为何将数据存储为栅格呢？有时只能将数据存储为栅格；例如，影像仅以栅格形式提供。然而，许多其他要素（例如点要素）和测量值（例如降雨量）既可以存储为栅格数据类型也可以存储为要素（矢量）数据类型。

GDAL

GDAL是栅格和矢量地理空间数据格式的转换器库，由开源地理空间基金会根据X / MIT样式的开源许可证发布。作为一个库，它为调用的应用程序提供了所有支持格式的单个栅格抽象数据模型和单个矢量抽象数据模型。它还带有用于数据转换和处理的各种有用的命令行实用程序。

栅格矢量化的流程：

1. 载入栅格

   首先加载需要矢量化的栅格数据集，读取栅格数据集的波段信息和空间信息。
   • 波段信息是矢量化的依据
   • 读取空间信息则是为了保证创建的矢量数据集和栅格数据集的坐标系保持一致

   栅格数据集如果没有空间信息，则矢量化的结果会出现上下反转的现象。

   Dataset dataset = gdal.Open(inRaster, gdalconstConstants.GA_ReadOnly);
   Band band = dataset.GetRasterBand(1);//栅格转矢量需要的波段信息
   SpatialReference prj = new SpatialReference();
   if (!dataset.GetProjectionRef().isEmpty()) {
       prj.ImportFromWkt(dataset.GetProjectionRef());//栅格数据的坐标系作为矢量化后的坐标系
   }
   

2. 创建矢量

   创建一个新的、空的矢量数据集来接收栅格矢量化后的数据

   String fileName = outShp.substring(outShp.lastIndexOf("\\") + 1);//带后缀的文件名
   String name = fileName.substring(0, fileName.lastIndexOf("."));//不带后缀
   Layer layer = ShpProcessing.createLayer("ESRI Shapefile", outShp, name, prj);
   FieldDefn field = new FieldDefn("value", ogr.OFTInteger);//创建一个字段用来存储栅格的像素值
   layer.CreateField(field, 1);
   

3. 矢量化

   执行gdal内置的的矢量化方法：gdal.Polygonize()

   gdal.Polygonize(band, null, layer, 0);
   layer.SyncToDisk();
   

根据指定坐标点读取栅格数据像元值：

• 根据文件路径获取到对应的DataSet类
• 从DataSet类获取Band类
• 调用Band类的ReadRaster获取一个一维数组(这个一维数组就是该栅格数据的某一区域的值。当ReadRaster的参数范围足够小时，就可以获取到一个像元的值，即一维数组只有一个值)

由于ReadRaster的参数需要像素坐标来定位，所以需要先使用DataSet类的GetGeoTransform方法获取像素坐标转地图坐标的参数，然后根据像素坐标转地图坐标的公式和指定的地图坐标点，解二元一次方程组，得到像素坐标。从而调用ReadRaster方法

package com.myself.raster;
 
import org.gdal.gdal.Band;
import org.gdal.gdal.Dataset;
import org.gdal.gdal.Driver;
import org.gdal.gdal.gdal;
import org.gdal.gdalconst.gdalconstConstants;
 
import java.util.Arrays;
 
public class RasterDemo {
 
    public static void main(String[] args) {
        String filepath = "F:\\raster\\Demo.tif";
        //注册
        gdal.AllRegister();
        //打开文件获取数据集
        Dataset dataset = gdal.Open(filepath,
                gdalconstConstants.GA_ReadOnly);
        if (dataset == null) {
            System.out.println("打开"+filepath+"失败"+gdal.GetLastErrorMsg());
            System.exit(1);
        }
        //获取驱动
        Driver driver = dataset.GetDriver();
        //获取驱动信息
        System.out.println("driver long name: " + driver.getLongName());
        //获取栅格数量
        int bandCount = dataset.getRasterCount();
        System.out.println("RasterCount: " + bandCount);
        //构造仿射变换参数数组，并获取数据
        double[] gt = new double[6];
        dataset.GetGeoTransform(gt);
        System.out.println("仿射变换参数"+Arrays.toString(gt));
 
        //指定经纬度
        double Latitude  = 86.053;
        double longitude = 16.529;
 
        //经纬度转换为栅格像素坐标
        int[] ColRow=Coordinates2ColRow(gt,Latitude,longitude);
 
        //判断是否坐标超出范围
        if(ColRow[0]<0||ColRow[1]<0||ColRow[0]>dataset.getRasterXSize()||ColRow[1]>dataset.getRasterYSize()){
            System.out.println(Arrays.toString(ColRow)+"坐标值超出栅格范围！");
            return;
        }
 
        //遍历波段，获取该点对应的每个波段的值并打印到屏幕
        for (int i = 0;i<bandCount;i++){
            Band band = dataset.GetRasterBand(i+1);
            double[] values = new double[1];
            band.ReadRaster(ColRow[0], ColRow[1], 1, 1, values);
            double value = values[0];
            System.out.println("横坐标："+Latitude +","+"纵坐标:"+longitude);
            System.out.format("Band"+(i+1)+": %s", value);
        }
        //释放资源
        dataset.delete();
 
    }
 
    /**
     * 将地图坐标转换为栅格像素坐标
     * @param gt 仿射变换参数
     * @param X 横坐标
     * @param Y 纵坐标
     * @return
     */
    public static int[] Coordinates2ColRow(double[] gt, double X, double Y){
        int[] ints = new int[2];
//        X = gt[0] + Xpixel*gt[1] + Yline*gt[2];
//        Y = gt[3] + Xpixel*gt[4] + Yline*gt[5];
//        消元法解二元一次方程组
//        X-gt[0]=Xpixel*gt[1] + Yline*gt[2]
//        Xpixel = (X-gt[0] - Yline*gt[2])/gt[1]
//        Y - gt[3] = ((X-gt[0] - Yline*gt[2])/gt[1])gt[4] + Yline*gt[5]
//        (Y - gt[3])*gt[1] = ((X-gt[0] - Yline*gt[2]))*gt[4] + Yline*gt[5]*gt[1]
//        (Y - gt[3])*gt[1] =(X-gt[0])*gt[4] - Yline*gt[2]*gt[4] + Yline*gt[5]*gt[1]
//        (Y - gt[3])*gt[1] - (X-gt[0])*gt[4] = Yline(gt[5]*gt[1]-gt[2]*gt[4])
 
        //向下取整,如果向上取整会导致计算结果偏大，从而在后面读取到邻近像元的数据
        double  Yline = Math.floor(((Y - gt[3])*gt[1] - (X-gt[0])*gt[4]) / (gt[5]*gt[1]-gt[2]*gt[4]));
        double  Xpixel = Math.floor((X-gt[0] - Yline*gt[2])/gt[1]);
        ints[0] = new Double(Xpixel).intValue();
        ints[1] = new Double(Yline).intValue();
        return ints;
    }
 
}

---

示例：

import java.io.File;
import org.gdal.gdal.Band;
import org.gdal.gdal.Dataset;
import org.gdal.gdal.Driver;
import org.gdal.gdal.gdal;
import org.gdal.gdalconst.gdalconstConstants;

public class gdaltest {
  public void dog() {
    gdal.AllRegister();
    String rasterFilePath ="F:\\倾斜摄影\\正射影像\\01可见光\\map\\result.tif";//测试用文件路径
    Dataset dataset = gdal.Open(rasterFilePath,
      gdalconstConstants.GA_ReadOnly);
    if (dataset == null) {
      System.out.println(gdal.GetLastErrorMsg());
    }
    Driver driver = dataset.GetDriver();
    System.out.println("driver short name: " + driver.getShortName());
    System.out.println("driver long name: " + driver.getLongName());
    System.out.println("metadata list: " + driver.GetMetadata_List());

    String proj = dataset.GetProjection();
    Band band = dataset.GetRasterBand(1);
    System.out.println(proj);
    System.out.println(band);
  }
  public static void main(String[] args) {
    gdaltest test = new gdaltest();
    test.dog();
  }
}

GeoTools

在地理信息系统的世界里有一类很重要的数据类型，那便是栅格数据。GeoTools对于栅格数据的支持并没有gdal强大，但既然他作为GeoTools的一部分，我们还是又必要了解一下它。

GeoTools支持的栅格数据有： arcgrid、geotiff、 grassraster、 image ( JPEG TIFF GIF PNG)、 imageio-ext-gdal、 imagemosaic、 imagepyramid、 JP2K、matlab

GeoTools支持的矢量数据有： wkt、csv、wkb、geojson、shapefile

GeoTools对于栅格数据的支持主要又由GridCoverage实现的。作为一名程序员，我们习惯于处理诸如JPEG、GIF或者PNG等格式的栅格数据。在地理空间方面，有一个Coverage个概念，他是空间定位要素的集合。非正式地，我们将地图和Coverage视为等同，当然，这是相对于地理意思而非编程思维上。

1. 栅格数据的读取（TIFF）

   public static Coverage readTiff(String tiffPath) throws IOException {
        File f = new File(tiffPath);
        ParameterValue<OverviewPolicy> policy = AbstractGridFormat.OVERVIEW_POLICY.createValue();
        policy.setValue(OverviewPolicy.IGNORE);
        ParameterValue<String> gridsize = AbstractGridFormat.SUGGESTED_TILE_SIZE.createValue();
        ParameterValue<Boolean> useJaiRead = AbstractGridFormat.USE_JAI_IMAGEREAD.createValue();
        useJaiRead.setValue(true);
        GridCoverage2D image = new GeoTiffReader(f).read(new GeneralParameterValue[]{policy, gridsize, useJaiRead});
        return image;
    }
   

2. 栅格数据的生成（TIFF）

   File file = new File(outTiffPath);
   GeoTiffWriter geoTiffWriter = new GeoTiffWriter(file);
   final GeoTiffFormat format = new GeoTiffFormat();
   final GeoTiffWriteParams wp = new GeoTiffWriteParams();
   // 设置写出参数
   wp.setCompressionMode(GeoTiffWriteParams.MODE_DEFAULT);
   wp.setTilingMode(GeoToolsWriteParams.MODE_DEFAULT);
   ParameterValueGroup paramWrite = format.getWriteParameters();
   paramWrite.parameter(AbstractGridFormat.GEOTOOLS_WRITE_PARAMS.getName().toString()).setValue(wp);
   geoTiffWriter.write((GridCoverage) coverage, paramWrite.values().toArray(new GeneralParameterValue[4]) );
   geoTiffWriter.dispose();
   
   

3. 栅格数据的加载

   AbstractGridFormat format = GridFormatFinder.findFormat(rasterFile);
   // this is a bit hacky but does make more geotiffs work
   Hints hints = new Hints();
   if (format instanceof GeoTiffFormat) {
       hints = new Hints(Hints.FORCE_LONGITUDE_FIRST_AXIS_ORDER, Boolean.TRUE);
   }
   reader = format.getReader(rasterFile, hints);
   

4. 波段提取

   GridCoverage2D cov = reader.read(null);
   int numBands = cov.getNumSampleDimensions();
   for (int i = 0; i < numBands; i++) {
       GridSampleDimension dim = cov.getSampleDimension(i);
       //这里依次输出 RED_BAND  GREEN_BANK BLUE_BANK
       System.out.println(dim.getDescription().toString());
   }
   

5. 样式创建

6. 读取tiff

   Java中如果要解析tiff，其实很多时候，我们都选择gdal的Java库来实现，毫无疑问，gdal确实在GIS数据处理方面非常的强悍，其实Geotools中很多有关栅格数据的解析，也是基于Gdal进行的封装，今天就简单了解使用GeoTools来解析Tiff数据。

   package com.dudu.gis;
   
   import org.geotools.gce.geotiff.GeoTiffReader;
   import org.geotools.geometry.GeneralEnvelope;
   import org.junit.Test;
   import org.opengis.coverage.grid.GridEnvelope;
   import org.opengis.referencing.crs.CoordinateReferenceSystem;
   
   import javax.media.jai.ImageLayout;
   import java.awt.color.ColorSpace;
   import java.awt.image.ColorModel;
   import java.awt.image.SampleModel;
   import java.io.IOException;
   
   /**
    * 一般对于tiff数据的读取，都会借助于gdal
    * Java API库，但其实Java的geotools也
    * 封装了读取tiff文件的API
    */
   public class GeoTiffReaderTest {
   
       private static final String TIF_PATH = "D:\\Gis开发\\数据\\影像数据\\tiff\\china.tif";
   
       @Test
       public void readTiffFile() throws IOException {
           /**
            * 使用GeoTiffReader
            */
           GeoTiffReader geoTiffReader = new GeoTiffReader(TIF_PATH);
           //获取到其bbox
           GeneralEnvelope envelope = geoTiffReader.getOriginalEnvelope();
           //获取到投影信息
           CoordinateReferenceSystem crs = geoTiffReader.getCoordinateReferenceSystem();
           //获取tiff的gridrange(网格范围，按像素说就是像素数量，也就是是tiff的rectangle)
           GridEnvelope gridEnvelope = geoTiffReader.getOriginalGridRange();
           //获取其图像相关信息
           ImageLayout imageLayout = geoTiffReader.getImageLayout();
           //获取colormodel
           ColorModel colorModel = imageLayout.getColorModel(null);
           //查看其通道数量
           int num = colorModel.getNumComponents();
           //是否支持透明
           boolean isAlpha = colorModel.hasAlpha();
           //获取颜色空间
           ColorSpace clorspace = colorModel.getColorSpace();
           //获取samplemodel
           SampleModel sampleModel = imageLayout.getSampleModel(null);
           //数据类型
           int datatype = sampleModel.getDataType();
       }
   }
   

QGIS

QGIS下载地址：https://www.qgis.org/zh-Hans/site/forusers/download.html