1 地理空间矢量实体的标准表示
地理几何实体（点、线、面、多点、多线、复杂多边形等）的规范表示是不同源数据交互的基础，因地理空间实体的复杂性，在实际业务中，空间实体格式的可视化解析，或计算异常是一件高频发生的事情，其多数情况下却只是格式不规范/不统一造成的。统一规范且合法地表示地理空间实体是必要的。
1，地理几何实体通用数据格式
1.1 WKT与GeoJson两种常用地理空间实体的简介
（1）WKT格式--OGC的标准
WKT(Well-known text)是开放地理空间联盟OGC（Open GIS Consortium ）制定的一种文本标记语言，用于表示矢量几何对象、空间参照系统及空间参照系统之间的转换。
（2）GeoJson格式--同时表示空间数据和属性数据的集合
GeoJSON 一种JSON格式的Feature信息输出格式，它便于被JavaScript等脚本语言处理，OpenLayers等地理库便是采用GeoJSON格式。
WKT与geojson可以表示的几何对象是一致的，包括点、线、面、几何集合四种：
1.Point, MultiPoint
2.LineString, MultiLineString
3.Polygon, MultiPolygon
4.GeometryCollection
1.2 WKT与GeoJson格式互转
/**
* wkt和geoJson互相转的工具类
*/
public class App {
private static WKTReader reader = new WKTReader();
private static final String GEO_JSON_TYPE="GeometryCollection";
private static final String WKT_TYPE="GEOMETRYCOLLECTION";
public static void main(String[] args) {
String wkt = "GEOMETRYCOLLECTION(POINT(4 6),LINESTRING(4 6,7 10))";
String wkt0="POLYGON((1 1,5 1,5 5,1 5,1 1),(2 2,2 3,3 3,3 2,2 2))";
JSONObject jsonObject = wktToJson(wkt);
System.out.println(jsonObject);
String s = jsonToWkt(jsonObject);
System.out.println("s = " + s);
}
/**
* wkt转Json
* @param wkt
* @return
*/
public static JSONObject wktToJson(String wkt) {
String json = null;
JSONObject jsonObject=new JSONObject();
try {
Geometry geometry = reader.read(wkt);
StringWriter writer = new StringWriter();
GeometryJSON geometryJSON=new GeometryJSON();
geometryJSON.write(geometry,writer);
json = writer.toString();
jsonObject = JSONObject.parseObject(json);
} catch (Exception e) {
System.out.println("WKT转GeoJson出现异常");
e.printStackTrace();
}
return jsonObject;
}
/**
* geoJson转wkt
* @param jsonObject
* @return
*/
public static String jsonToWkt(JSONObject jsonObject) {
String wkt = null;
String type = jsonObject.getString("type");
GeometryJSON gJson = new GeometryJSON();
try {
// {"geometries":[{"coordinates":[4,6],"type":"Point"},{"coordinates":[[4,6],[7,10]],"type":"LineString"}],"type":"GeometryCollection"}
if(App.GEO_JSON_TYPE.equals(type)){
// 由于解析上面的json语句会出现这个geometries属性没有采用以下办法
JSONArray geometriesArray = jsonObject.getJSONArray("geometries");
// 定义一个数组装图形对象
int size = geometriesArray.size();
Geometry[] geometries=new Geometry[size];
for (int i=0;i<size;i++){
String str = geometriesArray.get(i).toString();
// 使用GeoUtil去读取str
Reader reader = GeoJSONUtil.toReader(str);
Geometry geometry = gJson.read(reader);
geometries[i]=geometry;
}
GeometryCollection geometryCollection = new GeometryCollection(geometries,new GeometryFactory());
wkt=geometryCollection.toText();
}else {
Reader reader = GeoJSONUtil.toReader(jsonObject.toString());
Geometry read = gJson.read(reader);
wkt=read.toText();
}
} catch (IOException e){
System.out.println("GeoJson转WKT出现异常");
e.printStackTrace();
}
return wkt;
}
}
<geotools.version>17.1</geotools.version>
<dependency>
<groupId>org.geotools</groupId>
<artifactId>gt-geojson</artifactId>
<version>${geotools.version}</version>
</dependency>
<dependency>
<groupId>org.geotools</groupId>
<artifactId>gt-geometry</artifactId>
<version>${geotools.version}</version>
</dependency>
<dependency>
<groupId>com.vividsolutions</groupId>
<artifactId>jts</artifactId>
<version>1.13</version>
</dependency>
<dependency>
<groupId>org.geotools</groupId>
<artifactId>gt-epsg-hsql</artifactId>
<version>${geotools.version}</version>
</dependency>
<dependency>
<groupId>org.geotools</groupId>
<artifactId>gt-jts-wrapper</artifactId>
<version>${geotools.version}</version>
</dependency>
<!-- https://mvnrepository.com/artifact/org.geotools/gt-main -->
<dependency>
<groupId>org.geotools</groupId>
<artifactId>gt-main</artifactId>
<version>${geotools.version}</version>
</dependency>
<dependency>
<groupId>org.geotools</groupId>
<artifactId>gt-opengis</artifactId>
<version>${geotools.version}</version>
</dependency>
<!-- 引入json处理包 -->
<dependency>
<groupId>com.alibaba</groupId>
<artifactId>fastjson</artifactId>
<version>1.2.47</version>
</dependency>
<dependency>
<groupId>org.geotools</groupId>
<artifactId>gt-referencing</artifactId>
<version>${geotools.version}</version>
</dependency>1.3 WKT格式的多边形几何实体数据校准
WKT格式的格式化修复--可修复多边形未闭合、前缀后缀缺失重复、拓扑交叉造成不合法等多种常见格式异常问题。
package com.alibaba.dataworks.geo;
import java.util.regex.Matcher;
import java.util.regex.Pattern;
import java.util.Collection;
import java.util.Iterator;
import com.vividsolutions.jts.geom.Geometry;
import com.vividsolutions.jts.geom.GeometryFactory;
import com.vividsolutions.jts.geom.LineString;
import com.vividsolutions.jts.geom.LinearRing;
import com.vividsolutions.jts.geom.Polygon;
import com.vividsolutions.jts.io.WKTReader;
import com.vividsolutions.jts.operation.polygonize.Polygonizer;
/**
* Created by sheyu on 2018/4/7.
*/
public class polygonCleanAndFormat {
public GeometryFactory geometryFactory = new GeometryFactory();
public String evaluate(String polygonWkt,String IsReturnMultiPoygon,String isRepairHoles){
String resultWkt="";
if(polygonWkt==null||"".equals(polygonWkt)){
return polygonWkt;
}
try{
if(polygonWkt.contains("),(")||polygonWkt.contains("), (")){
resultWkt=polygonRepair(polygonWkt,IsReturnMultiPoygon,isRepairHoles);
}else{
polygonWkt=removeLetter(polygonWkt);
polygonWkt=removeSpChar(polygonWkt);
if("polygon".equals(polygonWkt.substring(0,7))||"POLYGON".equals(polygonWkt.substring(0,7))){
polygonWkt=polygonWkt.replace("POLYGON((", "").replace("polygon((", "").replace("))", "")
.replace("POLYGON ((", "").replace(", ", ",");
}
String coordinateArray[]=polygonWkt.split(",");
if(!coordinateArray[0].equals(coordinateArray[coordinateArray.length-1])){
polygonWkt+=","+coordinateArray[0];
}
polygonWkt="POLYGON (("+polygonWkt+"))";
resultWkt=polygonRepair(polygonWkt,IsReturnMultiPoygon,isRepairHoles);
if("parms error".equals(resultWkt)){
String secPoint=(Double.parseDouble(coordinateArray[1].split("\\s+")[0])+0.00001)+" "+coordinateArray[1].split("\\s+")[1];
polygonWkt=polygonWkt.replace(coordinateArray[1], secPoint);
resultWkt=polygonRepair(polygonWkt,IsReturnMultiPoygon,isRepairHoles);
}
}
return resultWkt.replace("POLYGON ((", "POLYGON((").replace(", ", ",");
}catch (Exception e) {
e.printStackTrace();
return "parms error";
}
}
@SuppressWarnings("unchecked")
public String polygonRepair(String polygonWkt,String isReturnMultiPoygon,String isRepairHoles) {
String resultWkt="";
try{
WKTReader reader = new WKTReader( geometryFactory );
Polygon polygon = (Polygon) reader.read(polygonWkt);
if(polygon.isValid()){
polygon.normalize();
if(isRepairHoles.equals("true")||isRepairHoles.equals("TRUE")){
LinearRing[] holes=null;
polygon=new Polygon((LinearRing) polygon.getExteriorRing(),holes, geometryFactory);
}
resultWkt=polygon.toText();
}else{
Polygonizer polygonizer = new Polygonizer();
addPolygon(polygon, polygonizer);
Geometry geom=toPolygonGeometry(polygonizer.getPolygons(), polygon.getFactory());
if(isReturnMultiPoygon.equals("true")||isReturnMultiPoygon.equals("TRUE")){
resultWkt=geom.toText();
}else{
polygon=(Polygon) geom.getGeometryN(0);
double polygonArea=polygon.getArea();
for(int i=0;i<geom.getNumGeometries();i++){
if(geom.getGeometryN(i).getArea()>polygonArea){
polygon=(Polygon) geom.getGeometryN(i);
polygonArea=polygon.getArea();
}
}
if(isRepairHoles.equals("true")||isRepairHoles.equals("TRUE")){
polygon=new Polygon((LinearRing) polygon.getExteriorRing(), null, geometryFactory);
}
resultWkt=polygon.toText().replace("POLYGON ((", "POLYGON((");
}
}
}catch (Exception e) {
return "parms error";
}
return resultWkt;
}
void addPolygon(Polygon polygon, Polygonizer polygonizer){
addLineString(polygon.getExteriorRing(), polygonizer);
for(int n = polygon.getNumInteriorRing(); n-- > 0;){
addLineString(polygon.getInteriorRingN(n), polygonizer);
}
}
void addLineString(LineString lineString, Polygonizer polygonizer){
if(lineString instanceof LinearRing){ // LinearRings are treated differently to line strings : we need a LineString NOT a LinearRing
lineString = lineString.getFactory().createLineString(lineString.getCoordinateSequence());
}
// unioning the linestring with the point makes any self intersections explicit.
com.vividsolutions.jts.geom.Point point = lineString.getFactory().createPoint(lineString.getCoordinateN(0));
Geometry toAdd = lineString.union(point);
//Add result to polygonizer
polygonizer.add(toAdd);
}
Geometry toPolygonGeometry(Collection<Polygon> polygons, GeometryFactory factory){
switch(polygons.size()){
case 0:
return null; // No valid polygons!
case 1:
return polygons.iterator().next(); // single polygon - no need to wrap
default:
//polygons may still overlap! Need to sym difference them
Iterator<Polygon> iter = polygons.iterator();
Geometry ret = iter.next();
while(iter.hasNext()){
ret = ret.symDifference(iter.next());
}
return ret;
}
}
public static String removeLetter(String value){
Pattern p = Pattern.compile("[a-zA-z]");
Matcher matcher = p.matcher(value);
String result = matcher.replaceAll("");
return result;
}
public static String removeSpChar(String value){
String regEx="[\n`~!@#$%^&*()+=|{}':;'\\[\\]<>/?~！@#￥%……&*（）——+|{}【】‘；：”“’。，、？]";
Pattern p = Pattern.compile(regEx);
Matcher matcher = p.matcher(value);
String result = matcher.replaceAll("");
return result;
}
public static void main(String[] args) {
polygonCleanAndFormat ddd=new polygonCleanAndFormat();
String polygonWkt="OLYGON((121.133616 28.818738,121.132307 28.815693,121.134871 28.815354,121.13559 28.818184,121.133616 28.818738)";
String resultWkt=ddd.evaluate(polygonWkt, "false", "true");
//polygonWkt=removeLetter(polygonWkt);
//polygonWkt=removeSpChar(polygonWkt);
System.out.println(resultWkt);
}
}2，基础地理几何体拓扑关系及计算
拓扑关系（ topological relation），指满足拓扑几何学原理的各空间数据间的相互关系。即用结点、弧段、多边形和岛所表示的实体之间的邻接、关联、包含和连通关系。如：点与点的邻接关系、点与面的包含关系、线与面的相离关系、面与面的重合关系等。
非拓扑关系计算包括两点之间的距离； 一个点指向另一个点的方向；弧段的长度；一个区域的周长；一个区域的面积。
2.1 地理空间几何实体拓扑关系
地理空间实体点、线、面两两拓扑关系包括：相离（disjoint）、相接(meet)、相交(overlap)、相等(equal)、包含且边界相交(cover-by)、包含且边界不交(contains)、包含于且边界相交(cover)、包含于且边界不交(inside)，如下图所示：
所谓地理空间实体之间的拓扑计算包括两个方面，一是拓扑关系的判断，二是拓扑交互结果的计算，当拓扑实体尤其是多边形面的相交时，会有多种情况--不同方式相交的面积的计算：
地理空间实体的拓扑计算的代码实现可以参考开源库JTS（Java，对应c++版本为Geos），部分拓扑计算开源代码：
import org.geotools.geometry.jts.JTSFactoryFinder;
import com.vividsolutions.jts.geom.Coordinate;
import com.vividsolutions.jts.geom.Geometry;
import com.vividsolutions.jts.geom.GeometryFactory;
import com.vividsolutions.jts.geom.LineString;
import com.vividsolutions.jts.geom.Point;
import com.vividsolutions.jts.io.ParseException;
import com.vividsolutions.jts.io.WKTReader;
/**
* Class GeometryRelated.java
* Description 二元比较集合。二元比较以两个几何对象作为参数，返回一个Boolean类型的值，
* 来指明这两个几何对象是否具有指定的空间关系。支持的空间关系包括：
* equals、disjoint、intersects, touches, crosses, within, contains, overlaps
*/
public class GeometryRelated {
private GeometryFactory geometryFactory = JTSFactoryFinder.getGeometryFactory( null );
public Point createPoint(String lon,String lat){
Coordinate coord = new Coordinate(Double.parseDouble(lon), Double.parseDouble(lat));
Point point = geometryFactory.createPoint( coord );
return point;
}
/**
* will return true as the two line strings define exactly the same shape.
* 两个几何对象是否是重叠的
* @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 ==geometry2; false
//check if two geometries are exactly equal; right down to the coordinate level.
// return geometry1.equalsExact(geometry2); false
return geometry1.equals(geometry2);//true
}
/**
* The geometries have no points in common
* 几何对象没有交点(相邻)
* @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);
}
/**
* The geometries have at least one point in common.
* 至少一个公共点(相交)
* @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);
}
/**
* @param args
* @throws ParseException
*/
public static void main(String[] args) throws ParseException {
GeometryRelated gr = new GeometryRelated();
System.out.println(gr.equalsGeo());
System.out.println(gr.disjointGeo());
System.out.println(gr.intersectsGeo());
}
}2.2 地理空间几何实体空间属性计算
如果说地理实体之间的拓扑计算为相互关系的计算，那么非拓扑计算则是独立实体的属性计算，如距离、面积、周长等，如果独立实体属性计算的坐标是经纬度，那么就要先进行投影到平面后再计算，参考代码：
package com.alibaba.dataworks.geo;
import org.geotools.geometry.jts.JTS;
import org.geotools.referencing.CRS;
import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.geom.GeometryFactory;
import org.locationtech.jts.geom.Polygon;
import org.locationtech.jts.geom.Point;
import org.locationtech.jts.io.ParseException;
import org.locationtech.jts.io.WKTReader;
import org.opengis.referencing.FactoryException;
import org.opengis.referencing.crs.CoordinateReferenceSystem;
import org.opengis.referencing.operation.MathTransform;
import org.opengis.referencing.operation.TransformException;
/**
* Created by sheyu on 2021/4/7.
*/
public class MercatorReference extends UDF {
public static GeometryFactory geometryFactory = new GeometryFactory();
//计算多边形投影面积、周长及投影后坐标本身
public String evaluate(String polygonWkt,String type,String srid) {
String reverseWkt=lntlatReverse(polygonWkt,null);
double value=0.0;
try{
Polygon geom=createPolygonByWKT(reverseWkt);
if(srid==null||"".equals(srid)){
double centerLine= Double.parseDouble(reverseWkt.split(",")[0].split(" ")[1]);//用于计算中央子午线参考经度
srid=String.valueOf(4488+Math.round(centerLine/3));
}
Geometry geometryMercator=transform(geom,srid);
if("area".equals(type)){
value=geometryMercator.getArea();
}else if("perimeter".equals(type)){
value=geometryMercator.getLength();
}else{
return geometryMercator.toText().replace(", ", ",").replace("POLYGON ((", "POLYGON((");
}
return String.valueOf(value);
}catch (Exception e) {
e.printStackTrace();
}
return "param error";
}
//计算单点的投影坐标
public String evaluate(String lng,String lat) {
String point="POINT("+lat+" "+lng+")";
try{
Point geom=createPointByWKT(point);
double centerLine= Double.parseDouble(lng);//用于计算中央子午线参考经度
String srid=String.valueOf(4488+Math.round(centerLine/3));
Geometry geometryMercator=transform(geom,srid);
return geometryMercator.toText().replace("POINT (", "POINT(");
}catch (Exception e) {
e.printStackTrace();
}
return "param error";
}
//计算两点球面距离--投影后距离
public String evaluate(String lng1,String lat1,String lng2,String lat2,String srid) {
String point1="POINT("+lat1+" "+lng1+")";
String point2="POINT("+lat2+" "+lng2+")";
try{
Point geom1=createPointByWKT(point1);
Point geom2=createPointByWKT(point2);
if(srid==null||"".equals(srid)){
double centerLine= (Double.parseDouble(lng1)+Double.parseDouble(lng2))/2;//用于计算中央子午线参考经度
srid=String.valueOf(4488+Math.round(centerLine/3));
}
Geometry geometryMercator1=transform(geom1,srid);
Geometry geometryMercator2=transform(geom2,srid);
return String.valueOf(geometryMercator1.distance(geometryMercator2));
}catch (Exception e) {
e.printStackTrace();
}
return "param error";
}
//计算距离--投影后距离
public String evaluate(String wkt1,String wkt2,String type1,String type2) {
String reversewkt1=lntlatReverse(wkt1,null);
String reversewkt2=lntlatReverse(wkt2,null);
try{
Geometry geom1=null;Geometry geom2=null;
double centerLine=120.0;
if("1".equals(type1)){
reversewkt1=lntlatReverse(wkt1,"point");
geom1=createPointByWKT(reversewkt1);
centerLine=Double.parseDouble(wkt1.split("\\(")[1].split(" ")[0]); //用于计算中央子午线参考经度
}else if("2".equals(type1)){
geom1=createPolygonByWKT(reversewkt1);
centerLine=Double.parseDouble(reversewkt1.split(",")[0].split(" ")[1]);
}
if("1".equals(type2)){
reversewkt2=lntlatReverse(wkt2,"point");
geom2=createPointByWKT(reversewkt2);
}else if("2".equals(type2)){
geom2=createPolygonByWKT(reversewkt2);
}
String srid=String.valueOf(4488+Math.round(centerLine/3));
Geometry geometryMercator1=transform(geom1,srid);
Geometry geometryMercator2=transform(geom2,srid);
return String.valueOf(geometryMercator1.distance(geometryMercator2));
}catch (Exception e) {
e.printStackTrace();
}
return "param error";
}
public static Geometry transform(Geometry geom,String srid) throws FactoryException,TransformException {
// WGS84(一般项目中常用的是CSR:84和EPSG:4326)
CoordinateReferenceSystem sourceCRS = CRS.decode("EPSG:4326");
// Pseudo-MercatorReference(墨卡托投影)
CoordinateReferenceSystem targetCRS = CRS.decode("EPSG:"+srid);
MathTransform transform = CRS.findMathTransform(sourceCRS, targetCRS, true);
Geometry geometryMercator = JTS.transform(geom, transform);
// 面积、周长
// System.out.println(geometryMercator.getArea());
// System.out.println(geom.getArea()*1e10);
// System.out.println(geometryMercator.getLength());
return geometryMercator;
}
public static Polygon createPolygonByWKT(String Polygon_points) throws ParseException{
WKTReader reader = new WKTReader( geometryFactory );
Polygon polygon = (Polygon) reader.read(Polygon_points);//"POLYGON((20 10, 30 0, 40 10, 30 20, 20 10))";
return polygon;
}
public static Point createPointByWKT(String points) throws ParseException{
WKTReader reader = new WKTReader( geometryFactory );
Point point = (Point) reader.read(points);//"POLYGON((20 10, 30 0, 40 10, 30 20, 20 10))";
return point;
}
private static String lntlatReverse(String polygonWkt,String type) {
if(polygonWkt==null||"".equals(polygonWkt)){
return polygonWkt;
}
String resultWkt="POLYGON((";
String polygonWktArray[]=polygonWkt.replace("POINT (", "").replace("POINT(", "").replace("POLYGON((", "").replace("POLYGON ((", "").replace(")", "").replace(" ", ",").split(",");
for(int i=0;i<polygonWktArray.length;i+=2){
if(i<polygonWktArray.length-2){
resultWkt+=polygonWktArray[i+1]+" "+polygonWktArray[i]+",";
}else{
resultWkt+=polygonWktArray[i+1]+" "+polygonWktArray[i]+"))";
}
}
if("point".equals(type)){
resultWkt=resultWkt.replace("POLYGON(","POINT").replace("))", ")");
}
return resultWkt;
}
public static void main(String[] args) {
String polygonWkts="POLYGON((116.700369 39.573407,116.701743 39.541908,116.746031 39.546144,116.73659 39.575656,116.700369 39.573407))";
String point="POINT(116.700369 39.573407)";
MercatorReference reference=new MercatorReference();
try{
String value=reference.evaluate(polygonWkts, "area", null);
System.out.println(value);
value=reference.evaluate(polygonWkts, "perimeter", null);
System.out.println(value);
value=reference.evaluate(polygonWkts, null, null);
System.out.println(value);
value=reference.evaluate(point,polygonWkts, "1","2");
System.out.println(value);
}catch (Exception e) {
e.printStackTrace();
}
}
}其他更多非拓扑计算参考个人先前ATA文章：loading
2，场景化通用高阶几何实体拓扑计算
在业务应用中，经常也会有基于场景的地理几何实体之间的复杂计算，如快递外卖场景中的快递员配送区分配、任务拆分、算法场景对地理空间的细粒度格网表征需求等。其需要对地理空间多实体（多边形）的无缝融合、智能拆分，geohash格网化计算等。
2.1 AOI（多边形）的无缝融合
AOI（多边形）的无缝融合算法是polygon buffer和union算法的融合算法，具体实现如下：
package com.cainiao.cnlbs.udf.odpsudf.sheyu;
import com.aliyun.odps.udf.UDF;
import com.vividsolutions.jts.geom.Geometry;
import com.vividsolutions.jts.geom.GeometryCollection;
import com.vividsolutions.jts.geom.GeometryFactory;
import com.vividsolutions.jts.geom.Polygon;
import com.vividsolutions.jts.io.ParseException;
import com.vividsolutions.jts.io.WKTReader;
import com.vividsolutions.jts.operation.buffer.BufferOp;
import com.vividsolutions.jts.operation.buffer.BufferParameters;
import com.cainiao.lbs.geo.PartitionTopologyOptimization;
/**
* Created by sheyu on 2018/8/29.
*/
public class PolygonGroupUnion extends UDF {
public static GeometryFactory geometryFactory = new GeometryFactory();
//polygonWkts：wkt格式，分号分隔;tolerateDistance：融合距离;isRepairHoles：是否修复空洞
public String evaluate(String polygonWkts,String tolerateDistance,String isRepairHoles){
String resultWkt="";
try{
if(polygonWkts==null||"".equals(polygonWkts)){
return resultWkt;
}
String polygonWktArray[]=polygonWkts.split(";");
Geometry[] revGeoms = new Geometry[polygonWktArray.length];
for(int i=0;i<polygonWktArray.length;i++){
revGeoms[i]=createPolygonByWKT(polygonWktArray[i]);
}
GeometryCollection collection=geometryFactory.createGeometryCollection(revGeoms);
double bufferDistance=tolerateDistance==null?0.0003:Double.parseDouble(tolerateDistance);
BufferParameters param=new BufferParameters(16,3,2,5);
Geometry unionGeom= BufferOp.bufferOp(BufferOp.bufferOp(collection,bufferDistance,param),-1*bufferDistance,param).union();
resultWkt=PartitionTopologyOptimization.getMaxAreaPolygon(unionGeom).toText();
if("true".equals(isRepairHoles)||"1".equals(isRepairHoles)){
polygonCleanAndFormat polygonClean=new polygonCleanAndFormat();
resultWkt=polygonClean.evaluate(resultWkt, "true", "true");
}
return resultWkt;
}catch (Exception e) {
e.printStackTrace();
return "parms error";
}
}
public static Polygon createPolygonByWKT(String Polygon_points) throws ParseException{
WKTReader reader = new WKTReader( geometryFactory );
Polygon polygon = (Polygon) reader.read(Polygon_points);//"POLYGON((20 10, 30 0, 40 10, 30 20, 20 10))";
return polygon;
}
}2.2 AOI（多边形）的Geohash格网拆分
AOI（多边形）的格网拆分算法，具体实现如下：
package com.cainiao.lbs.geo;
import ch.hsr.geohash.BoundingBox;
import java.util.ArrayList;
import java.util.List;
import ch.hsr.geohash.GeoHash;
import com.vividsolutions.jts.geom.Coordinate;
import com.vividsolutions.jts.geom.Geometry;
import com.vividsolutions.jts.geom.GeometryFactory;
import com.vividsolutions.jts.geom.Polygon;
import com.vividsolutions.jts.io.WKTReader;
import com.vividsolutions.jts.io.WKTWriter;
import com.aliyun.odps.udf.UDF;
/**
* Created by sheyu on 2018/2/27.
*/
public class PolygonToGeoHash extends UDF{
List<String> resultList = new ArrayList<>();
private GeometryFactory geometryFactory = new GeometryFactory();
public String evaluate(String wkt,String lentghOfGeoHash) {
String result="";
if(resultList!=null){
resultList.clear();
}
try{
WKTReader reader = new WKTReader(geometryFactory);
Polygon polygon = (Polygon) reader.read(wkt);
int lentgh=Integer.parseInt(lentghOfGeoHash);
List<String> resultGeohashList=getGeohahList(polygon,lentgh);
if(resultGeohashList!=null){
for(int i=0;i<resultGeohashList.size();i++){
result+=resultGeohashList.get(i);
if(i<resultGeohashList.size()-1){
result+=",";
}
}
}
return result;
}catch (Exception e) {
e.printStackTrace();
return "parms error";
}
}
public List<String> getGeohahList(Polygon polygon, int lentghOfGeoHash){
if( polygon==null|| lentghOfGeoHash<0 || lentghOfGeoHash>12){
resultList.clear();
return resultList;
}
BoundingBox boundingBox = getPolygonBoundingBox( polygon);
GeoHash hashUpperLeft = GeoHash.withCharacterPrecision( boundingBox.getUpperLeft().getLatitude(),
boundingBox.getUpperLeft().getLongitude(),lentghOfGeoHash);
GeoHash hashLowerRight = GeoHash.withCharacterPrecision( boundingBox.getLowerRight().getLatitude(),
boundingBox.getLowerRight().getLongitude(),lentghOfGeoHash);
if (hashLowerRight.equals( hashUpperLeft )){
resultList.add( hashLowerRight.toBase32());
System.out.println( transeGeohashToWKT(hashLowerRight) +";");
return resultList;
}
double perLng = hashUpperLeft.getBoundingBox().getLongitudeSize();
double perLat = hashUpperLeft.getBoundingBox().getLatitudeSize();
int lngSteps = (int) Math.ceil( ( boundingBox.getMaxLon() - boundingBox.getMinLon()) / perLng)+1;
int latSteps = (int) Math.ceil( ( boundingBox.getMaxLat() - boundingBox.getMinLat()) / perLat)+1 ;
if((lngSteps * latSteps) >10000){
// System.out.println("too much geohash.");
return resultList;
}
GeoHash lowerGeohash = hashUpperLeft;
for(int y=0; y<latSteps; y++){
GeoHash lineFirstHash = lowerGeohash;
for( int x=0; x<lngSteps; x++){
Polygon geohashpolygon=transeGeohashToPolygon(lineFirstHash);
if(polygon.intersects(geohashpolygon)==true){
resultList.add( lineFirstHash.toBase32());
}
//System.out.println( transeGeohashToWKT(lineFirstHash) +";");
lineFirstHash = lineFirstHash.getEasternNeighbour();
}
lowerGeohash = lowerGeohash.getSouthernNeighbour();
}
return resultList;
}
public String transeGeohashToWKT( GeoHash geoHash ){
if( geoHash==null){
return "";
}
BoundingBox boundingBox= geoHash.getBoundingBox();
// Envelope envelope = new Envelope( boundingBox.getMinLon(), boundingBox.getMaxLon(),
// boundingBox.getMinLat(), boundingBox.getMaxLat());
Coordinate[] coordinates = new Coordinate[]{
new Coordinate( boundingBox.getMaxLon(), boundingBox.getMaxLat()),
new Coordinate( boundingBox.getMinLon(), boundingBox.getMaxLat()),
new Coordinate( boundingBox.getMinLon(), boundingBox.getMinLat()),
new Coordinate( boundingBox.getMaxLon(), boundingBox.getMinLat()),
new Coordinate( boundingBox.getMaxLon(), boundingBox.getMaxLat())
};
Geometry geometry = geometryFactory.createLineString(coordinates);
WKTWriter wktWriter = new WKTWriter();
return wktWriter.write( geometry);
}
public Polygon transeGeohashToPolygon( GeoHash geoHash ){
if( geoHash==null){
return null;
}
BoundingBox boundingBox= geoHash.getBoundingBox();
Coordinate[] coordinates = new Coordinate[]{
new Coordinate( boundingBox.getMaxLon(), boundingBox.getMaxLat()),
new Coordinate( boundingBox.getMinLon(), boundingBox.getMaxLat()),
new Coordinate( boundingBox.getMinLon(), boundingBox.getMinLat()),
new Coordinate( boundingBox.getMaxLon(), boundingBox.getMinLat()),
new Coordinate( boundingBox.getMaxLon(), boundingBox.getMaxLat())
};
Polygon polygon = geometryFactory.createPolygon(coordinates);
return polygon;
}
public BoundingBox getPolygonBoundingBox( Polygon polygon){
double minLng = Double.MAX_VALUE;
double maxLng = Double.MIN_VALUE;
double minLat = Double.MAX_VALUE;
double maxLat = Double.MIN_VALUE;
Coordinate[] coordinates = polygon.getCoordinates();
for( int i=0; i<coordinates.length; i++){
minLng = Math.min(minLng, coordinates[i].x);
maxLng = Math.max(maxLng, coordinates[i].x);
minLat = Math.min(minLat, coordinates[i].y);
maxLat = Math.max( maxLat, coordinates[i].y);
}
BoundingBox boundingBox = new BoundingBox(minLat, maxLat, minLng, maxLng);
return boundingBox;
}
}2.3 复杂AOI（多边形）的边界平滑
package com.cainiao.lbs.geo;
import com.aliyun.odps.udf.UDF;
import com.vividsolutions.jts.geom.*;
import com.vividsolutions.jts.io.WKTReader;
import com.vividsolutions.jts.operation.polygonize.Polygonizer;
import java.util.Collection;
import java.util.Iterator;
/**
* Created by sheyu on 2018/4/7.
*/
public class smoothPolygon extends UDF {
public GeometryFactory geometryFactory = new GeometryFactory();
public String evaluate(String polygonWkt,String angleLimit,String distanceLimit){
try{
Polygon polygon=CommonUtils.createPolygonByWKT(polygonWkt);
Polygon resultPolygon=IsConcave.smoothPolygon(polygon,Double.parseDouble(angleLimit),Double.parseDouble(distanceLimit));
return resultPolygon.toText().replace(", ", ",").replace("POLYGON ((", "POLYGON((");
}catch (Exception e) {
e.printStackTrace();
return "parms error";
}
}
}
//IsConcave.smoothPolygon方法
/**
* 多边形锯齿优化
* @param angleLimit 顶点角度阈值
* @param distanceLimit 边长距离阈值
* @param polygon 待优化polygon
* @return 优化后多边形
*/
public static Polygon smoothPolygon(Polygon polygon,double angleLimit,double distanceLimit) {
if(polygon==null||polygon.getNumPoints()<=3){
return polygon;
}
Coordinate[] cs = polygon.getCoordinates();
Coordinate[] polygonArray = polygon.getCoordinates();
if(distanceLimit==0){
distanceLimit=polygon.getLength()/polygon.getNumPoints()*0.05;
}else{
distanceLimit=polygon.getLength()/polygon.getNumPoints()*distanceLimit;
}
while(FirstPointCheck(polygonArray,distanceLimit,angleLimit)==false){
polygonArray=new Coordinate[polygonArray.length-1];
for(int i=0;i<polygonArray.length;i++){
if(i<polygonArray.length-1) {
polygonArray[i] = cs[i + 1];
}else{
polygonArray[i] = cs[1];
}
}
cs=polygonArray;
}
if(polygonArray.length<=4){
return new GeometryFactory().createPolygon(polygonArray);
}
int firstpointindex=0;
int lastpointindex=0;
int k=0;
for (int i = 0; i < cs.length-1; i++) {
lastpointindex=i+1;
if(cs[i].distance(cs[lastpointindex])<distanceLimit){
polygonArray[i]=null;
k++;
continue;
}
firstpointindex=i-1;
for(int h=firstpointindex;;h--){
if(h<0){
h=cs.length-2;
}
if(polygonArray[h]!=null){
firstpointindex=h;
break;
}
}
double angle=Math.abs(azimuthAngle(cs[i].x,cs[i].y, cs[firstpointindex].x,cs[firstpointindex].y)-azimuthAngle(cs[i].x,cs[i].y, cs[lastpointindex].x,cs[lastpointindex].y));
if (angle<angleLimit || 360-angle<angleLimit ){
polygonArray[i]=null;
k++;
}
}
Coordinate[] clist = new Coordinate[cs.length-k];
int h=0;
for(int i=0;i<polygonArray.length;i++){
if(polygonArray[i]!=null){
clist[h++]=polygonArray[i];
}
}
if(polygonArray[0]==null){
clist[h-1]=clist[0];
}
try{
return new GeometryFactory().createPolygon(clist);
}catch (Exception e) {
return polygon;
}
}2.4 其他更多复杂空间实体计算
其他更多实体计算，如相对geohash拆分更灵活的任意格网拆分、AOI（多边形）的聚类拆分、凹多边形的凸拆分等，参考个人之前ATA文章：loading
扫一扫
3259
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
