cutfilldemo/src/main/java/org/example/test/CreateRingGridFixGeo.java

864 lines
43 KiB
Java
Raw Normal View History

package org.example.test;
import org.example.alpha.AlphaShape;
import org.example.alpha.Vector2D;
import org.gdal.gdal.Dataset;
import org.gdal.gdal.Driver;
import org.gdal.gdal.WarpOptions;
import org.gdal.gdal.gdal;
import org.gdal.gdalconst.gdalconst;
import org.geotools.api.feature.simple.SimpleFeature;
import org.geotools.api.feature.simple.SimpleFeatureType;
import org.geotools.api.geometry.Bounds;
import org.geotools.api.geometry.Position;
import org.geotools.api.style.Style;
import org.geotools.coverage.grid.GridCoordinates2D;
import org.geotools.coverage.grid.GridCoverage2D;
import org.geotools.coverage.grid.GridGeometry2D;
import org.geotools.coverage.processing.Operations;
import org.geotools.feature.DefaultFeatureCollection;
import org.geotools.feature.simple.SimpleFeatureBuilder;
import org.geotools.feature.simple.SimpleFeatureTypeBuilder;
import org.geotools.gce.geotiff.GeoTiffReader;
import org.geotools.geojson.feature.FeatureJSON;
import org.geotools.geojson.geom.GeometryJSON;
import org.geotools.geometry.Position2D;
import org.geotools.geometry.jts.ReferencedEnvelope;
import org.geotools.map.FeatureLayer;
import org.geotools.map.Layer;
import org.geotools.map.MapContent;
import org.geotools.referencing.crs.DefaultGeographicCRS;
import org.geotools.styling.SLD;
import org.geotools.swing.JMapFrame;
import org.geotools.util.factory.Hints;
import org.locationtech.jts.geom.Point;
import org.locationtech.jts.geom.Polygon;
import org.locationtech.jts.geom.*;
import org.locationtech.jts.operation.buffer.BufferOp;
import org.locationtech.jts.operation.buffer.BufferParameters;
import org.springframework.beans.factory.annotation.Value;
import javax.media.jai.Interpolation;
import javax.media.jai.InterpolationBicubic2;
import java.awt.*;
import java.io.File;
import java.util.List;
import java.util.*;
public class CreateRingGridFixGeo {
/**
* 必须注册
*/
static {
gdal.AllRegister();
}
public static GridCoverage2D coverage;
public static String inputDemPath = "E:\\data\\CB测试\\geospatial_data_cloud\\ASTGTMV003_N30E119_dem.tif";
public static String outputDemPath = "E:\\data\\CB测试\\geospatial_data_cloud\\resample\\coverage_geo.tif";
public static double extendLen = 100;
public static GeometryFactory geometryFactory = new GeometryFactory();
public static List<Map<String,Object>> fillPoints;
public static List<Map<String,Object>> cutPoints;
public static double cutAllSum = 0;
public static double fillAllSum = 0;
public static double meter = 2; // 间距,单位为米。
public static double range = 10;//容差
public static boolean isProjcetion = false;
public static boolean isResample = true;
private static SimpleFeatureType createFeatureType() {
SimpleFeatureTypeBuilder builder = new SimpleFeatureTypeBuilder();
builder.setName("Location");
builder.setCRS(DefaultGeographicCRS.WGS84); // <- Coordinate reference system
// add attributes in order
builder.add("the_geom", Polygon.class);
builder.length(15).add("name", String.class); // <- 15 chars width for name field
builder.add("id", Integer.class);
// build the type
final SimpleFeatureType LOCATION = builder.buildFeatureType();
return LOCATION;
}
private static SimpleFeatureType createLineStringFeatureType() {
SimpleFeatureTypeBuilder builder = new SimpleFeatureTypeBuilder();
builder.setName("Location");
builder.setCRS(DefaultGeographicCRS.WGS84); // <- Coordinate reference system
// add attributes in order
builder.add("the_geom", LineString.class);
builder.length(15).add("name", String.class); // <- 15 chars width for name field
builder.add("id", Integer.class);
// build the type
final SimpleFeatureType LOCATION = builder.buildFeatureType();
return LOCATION;
}
private static SimpleFeatureType createPointFeatureType() {
SimpleFeatureTypeBuilder builder = new SimpleFeatureTypeBuilder();
builder.setName("Location");
builder.setCRS(DefaultGeographicCRS.WGS84); // <- Coordinate reference system
// add attributes in order
builder.add("the_geom", Point.class);
// builder.length(15).add("name", String.class); // <- 15 chars width for name field
builder.add("id", String.class);
builder.add("bottomHeight", double.class);
builder.add("topHeight", double.class);
builder.add("elev", double.class);
return builder.buildFeatureType();
}
public static void main(String[] args) throws Exception{
loadDem(inputDemPath);
GeometryFactory geometryFactory = new GeometryFactory();
List<double[]> coords = new ArrayList(){{
// 新给的dem,自己画的
// add(new double[]{469037.331529372138903, 3348374.219233838375658 });
// add(new double[]{469256.397423173999414, 3348360.942513002082705 });
// add(new double[]{469246.439882546663284, 3348171.749241081997752 });
// add(new double[]{469014.097267908335198, 3348194.983502545859665 });
// add(new double[]{469037.331529372138903, 3348374.219233838375658});
add(new double[]{119.678274321091934, 30.254540850428963 });
add(new double[]{119.680550940139568, 30.254426658580112 });
add(new double[]{119.680452999081325, 30.252719786115247 });
add(new double[]{119.678038175380664, 30.252923459684105 });
add(new double[]{119.678274321091934, 30.254540850428963 });
}};
// 创建矩形的四个角点坐标
Coordinate[] coordinates = new Coordinate[5];
coordinates[0] = new Coordinate(coords.get(0)[0], coords.get(0)[1]);
coordinates[1] = new Coordinate(coords.get(1)[0], coords.get(1)[1]);
coordinates[2] = new Coordinate(coords.get(2)[0], coords.get(2)[1]);
coordinates[3] = new Coordinate(coords.get(3)[0], coords.get(3)[1]);
// 闭合多边形,复制第一个坐标作为最后一个坐标
coordinates[4] = new Coordinate(coords.get(0)[0], coords.get(0)[1]);
// 计算设计标高,取四个角的高程平均值
// double centerHeight = getElevation(rectangle.getCentroid());
double corner1Height = getElevation( geometryFactory.createPoint(coordinates[0]));
double corner2Height = getElevation( geometryFactory.createPoint(coordinates[1]));
double corner3Height = getElevation( geometryFactory.createPoint(coordinates[2]));
double corner4Height = getElevation( geometryFactory.createPoint(coordinates[3]));
double targetHeight = (corner1Height+corner2Height+corner3Height+corner4Height)/4;
// double targetHeight = 158.5;
// double targetHeight = 145.7851371725257;
System.out.println("初始标高:"+targetHeight);
double startTime = System.currentTimeMillis();
Polygon rectangle = geometryFactory.createPolygon(coordinates);
// resampleDem(rectangle);
if(isResample){
// 重采样设置
BufferParameters bufferParameters = new BufferParameters();
bufferParameters.setEndCapStyle(BufferParameters.CAP_FLAT);
bufferParameters.setJoinStyle(BufferParameters.JOIN_MITRE);
BufferOp bufOp = new BufferOp(rectangle,bufferParameters);
Polygon bufferLinePolygon = (Polygon)bufOp.getResultGeometry(isProjcetion?extendLen:meter2degree(extendLen));
double resampleStartTime = System.currentTimeMillis();
double distance = isProjcetion?meter:meter2degree(meter);
reSample(bufferLinePolygon.getEnvelope(),inputDemPath,outputDemPath,distance,distance);
loadDem(outputDemPath);
System.out.println("重采样耗时:"+(System.currentTimeMillis()-resampleStartTime)+"毫秒");
}
double allSum = calcBlanceHeight(rectangle,targetHeight);
// double firstSum =0;
// double step = 2;//起始步长
// double nextHeight=targetHeight;
// int count = 1;
// while (Math.abs(allSum)-range>0){
// // if((firstSum>0&&allSum<0)||(firstSum<0&&allSum>0)){
// // fix: 采用线性方式去计算步长
// step = step*Math.abs(allSum)/Math.abs(firstSum-allSum);
// // }
// nextHeight = allSum<0?nextHeight-step:nextHeight+step;
// firstSum = calcBlanceHeight(rectangle,nextHeight);
// count++;
// System.out.println("firstSum:"+firstSum+"\t"+nextHeight);
// if(Math.abs(firstSum)-range<0){
// allSum = firstSum;
// break;
// }
// // if((firstSum>0&&allSum<0)||(firstSum<0&&allSum>0)){
// // fix: 采用线性方式去计算步长
// step = step*Math.abs(firstSum)/Math.abs(firstSum-allSum);
// // }
// double tmpHeight = firstSum<0?nextHeight-step:nextHeight+step;
// allSum = calcBlanceHeight(rectangle,tmpHeight);
// count++;
// if(Math.abs(tmpHeight-nextHeight)<0.0001){
// nextHeight = tmpHeight;
// break;
// }
// nextHeight = tmpHeight;
// System.out.println("allSum:"+allSum+"\t"+nextHeight);
// }
//
// // double radius = meter2degree(meter);
// double radius = isProjcetion?meter:meter2degree(meter);
// SimpleFeature cutBoundary = searchBourdary(cutPoints,radius,"挖方边界线");
// SimpleFeature fillBoundary = searchBourdary(fillPoints,radius,"填方边界线");
// FeatureJSON featureJSON = new FeatureJSON(new GeometryJSON(8));
//
// System.out.println("理想标高:"+nextHeight);
// System.out.println("总挖方量:"+cutAllSum);
// System.out.println("总填方量:"+fillAllSum);
// System.out.println("填挖方差值:"+allSum);
//
// System.out.println("单次耗时:"+(System.currentTimeMillis()-startTime)+"毫秒");
//
// System.out.println("总填方格子数量:"+fillPoints.size());
// System.out.println("总挖方格子数量:"+cutPoints.size());
//
// System.out.println("耗时:"+(System.currentTimeMillis()-startTime));
//
// System.out.println("挖方边界线:"+featureJSON.toString(cutBoundary));
// System.out.println("填方边界线:"+featureJSON.toString(fillBoundary));
// 输出:理想标高、总挖方量、总填方量、填挖方差值、边界数据
}
public static double calcBlanceHeight(Polygon rectangle ,double targetHeight) throws Exception{
// 创建GeometryFactory
List<Map<String,Object>> featuresProperties = new ArrayList<>();// 面要素属性集合
List<Map<String,Object>> tPointFeaturesProperties = new ArrayList<>();// 目标平面的点集属性集合
// 清空
fillPoints = new ArrayList<>();
cutPoints = new ArrayList<>();
cutAllSum = 0;
fillAllSum = 0;
// 间距,单位为米。
// double meter = 2;
double bottomRatio = 1.0/2;
double topRatio = 1.0/1.5;
// 创建多边形(矩形)
Map<String,Object> innerPolygon = new HashMap<>();
// Polygon rectangle = geometryFactory.createPolygon(coordinates);
innerPolygon.put("geom",rectangle);
innerPolygon.put("name","输入的面");
innerPolygon.put("id",1);
// BufferParameters bufferParameters = new BufferParameters();
// bufferParameters.setEndCapStyle(BufferParameters.CAP_FLAT);
// bufferParameters.setJoinStyle(BufferParameters.JOIN_MITRE);
// BufferOp bufOp = new BufferOp(rectangle,bufferParameters);
// Polygon bufferLinePolygon = (Polygon)bufOp.getResultGeometry(isProjcetion?extendLen:meter2degree(extendLen));
// innerPolygon.put("geom",bufferLinePolygon.getEnvelope());
featuresProperties.add(innerPolygon);
// 外接矩形打点
Map<String,Object> outPolygon1 = new HashMap<>();
Polygon out1 = (Polygon) rectangle.getEnvelope();//geometryFactory.createPolygon(cords4);
outPolygon1.put("geom",out1);
outPolygon1.put("name","输入面的外接矩形");
outPolygon1.put("id",6);
featuresProperties.add(outPolygon1);
Coordinate[] outCoords = out1.getCoordinates();
Coordinate[] outLine = new Coordinate[]{outCoords[0],outCoords[1]};
// 定义地图单位的距离
double distance = meter;
if(!isProjcetion){
distance = meter2degree(distance);
}
addPartPoint(outLine,geometryFactory,distance,rectangle,tPointFeaturesProperties,targetHeight,"0");
// 扩展面打点
List<List<Map<String,Object>>> pointRingFeaturesProperties = createRingPoint(rectangle,targetHeight,bottomRatio,topRatio,distance);
// 中心点
Point centroid = rectangle.getCentroid();
// 计算网格
Map<String,Object> result = calcGrid(pointRingFeaturesProperties,centroid,meter);
List<Map<String,Object>> fillPointList = (List<Map<String,Object>>)result.get("fillPointList");
List<Map<String,Object>> cutPointList = (List<Map<String,Object>>)result.get("cutPointList");
double fillSum = (double)result.get("fillSum");
double cutSum = (double)result.get("cutSum");
// 计算设计平面内的填挖方量
for(int i=0;i<tPointFeaturesProperties.size();i++){
Map<String, Object> pfp = tPointFeaturesProperties.get(i);
// double height = (double)pfp.get("height");
double height = (double)pfp.get("bottomHeight");
Point p = (Point)pfp.get("geom");
double elev = getElevation(p);
pfp.put("elev",elev);
if(elev>height){
cutPointList.add(pfp);
cutSum += (elev - height)*Math.pow(meter,2);
}else{
fillPointList.add(pfp);
fillSum += (elev - height)*Math.pow(meter,2);
}
}
System.out.println("总挖方量:"+cutSum);
System.out.println("总填方量:"+fillSum);
System.out.println("填挖方差量:"+(fillSum+cutSum));
// List<Map<String, Object>> boundary = searchBourdary(cutPointList, meter2degree(meter));
// Create a map content and add our shapefile to it
MapContent map = new MapContent();
map.setTitle("填挖方计算-创建网格");
// 构建添加面图层
final SimpleFeatureType TYPE =createFeatureType();
DefaultFeatureCollection featureCollection = new DefaultFeatureCollection();
SimpleFeatureBuilder featureBuilder = new SimpleFeatureBuilder(TYPE);
featuresProperties.forEach(props->{
featureBuilder.add(props.get("geom"));
featureBuilder.add(props.get("name"));
featureBuilder.add(props.get("id"));
SimpleFeature feature = featureBuilder.buildFeature((String.valueOf(props.get("id"))));
featureCollection.add(feature);
});
Style style = SLD.createSimpleStyle(featureCollection.getSchema());
Layer layer = new FeatureLayer(featureCollection, style);
map.addLayer(layer);
// 构建点要素图层
// List<List<Map<String,Object>>> pointList = new ArrayList(){{
//// add(result.get("fillPointList"));
// add(tPointFeaturesProperties);
// }};
pointRingFeaturesProperties.add(tPointFeaturesProperties);
List<List<Map<String,Object>>> pointList = pointRingFeaturesProperties;
pointList.stream().forEach(list->{
if(!list.isEmpty()){
final SimpleFeatureType PTYPE =createPointFeatureType();
DefaultFeatureCollection pointFeatureCollection = new DefaultFeatureCollection();
SimpleFeatureBuilder pointFeatureBuilder = new SimpleFeatureBuilder(PTYPE);
list.forEach(props->{
pointFeatureBuilder.add(props.get("geom"));
pointFeatureBuilder.add(props.get("id"));
pointFeatureBuilder.add(props.get("bottomHeight"));
pointFeatureBuilder.add(props.get("topHeight"));
pointFeatureBuilder.add(props.get("elev"));
SimpleFeature feature = pointFeatureBuilder.buildFeature((String.valueOf(props.get("id"))));
pointFeatureCollection.add(feature);
});
Style pointStyle = SLD.createSimpleStyle(pointFeatureCollection.getSchema());
Layer pointLlayer = new FeatureLayer(pointFeatureCollection, pointStyle);
map.addLayer(pointLlayer);
}
});
final SimpleFeatureType PTYPE =createPointFeatureType();
DefaultFeatureCollection pointFeatureCollection = new DefaultFeatureCollection();
SimpleFeatureBuilder pointFeatureBuilder = new SimpleFeatureBuilder(PTYPE);
List<Map<String,Object>> list = (List<Map<String,Object>>)result.get("cutPointList");
cutPointList.forEach(props->{
pointFeatureBuilder.add(props.get("geom"));
pointFeatureBuilder.add(props.get("id"));
pointFeatureBuilder.add(props.get("bottomHeight"));
pointFeatureBuilder.add(props.get("topHeight"));
SimpleFeature feature = pointFeatureBuilder.buildFeature((String.valueOf(props.get("id"))));
pointFeatureCollection.add(feature);
});
Style pointStyle = SLD.createSimpleStyle(pointFeatureCollection.getSchema(),Color.red);
Layer pointLlayer = new FeatureLayer(pointFeatureCollection, pointStyle);
map.addLayer(pointLlayer);
final SimpleFeatureType FillPTYPE =createPointFeatureType();
DefaultFeatureCollection pointFillFeatureCollection = new DefaultFeatureCollection();
SimpleFeatureBuilder pointFillFeatureBuilder = new SimpleFeatureBuilder(FillPTYPE);
// List<Map<String,Object>> fillList = (List<Map<String,Object>>)result.get("fillPointList");
fillPointList.forEach(props->{
pointFillFeatureBuilder.add(props.get("geom"));
pointFillFeatureBuilder.add(props.get("id"));
pointFillFeatureBuilder.add(props.get("bottomHeight"));
pointFillFeatureBuilder.add(props.get("topHeight"));
SimpleFeature feature = pointFillFeatureBuilder.buildFeature((String.valueOf(props.get("id"))));
pointFillFeatureCollection.add(feature);
});
Style pointFillStyle = SLD.createSimpleStyle(pointFillFeatureCollection.getSchema(),Color.blue);
Layer pointFillLlayer = new FeatureLayer(pointFillFeatureCollection, pointFillStyle);
map.addLayer(pointFillLlayer);
JMapFrame.showMap(map);
cutAllSum = cutSum;
cutPoints = cutPointList;
fillAllSum = fillSum;
fillPoints = fillPointList;
return fillSum + cutSum;
// return 0;
}
// private static List<Map<String, Object>> searchBourdary(List<Map<String, Object>> fillPointList, double radius,String name) {
//
// AlphaShape alphaShape = new AlphaShape(radius);
// List<Map<String, Object>> result = new ArrayList<>();
// // 构建顶点集合对象
// List<Vector2D> points = new ArrayList<>();
// fillPointList.forEach(pointMap->{
// Point p = (Point)pointMap.get("geom");
// Vector2D vector2D = new Vector2D(p.getX(), p.getY());
//
// points.add(vector2D);
// });
// List<List<Vector2D>> lists = alphaShape.boundaryPoints(points);
// lists.forEach(pointList->{
// pointList.forEach(point->{
// HashMap<String, Object> tmp = new HashMap<>();
// tmp.put("geom",geometryFactory.createPoint(new Coordinate(point.getX(),point.getY())));
// tmp.put("id",point.getIndex());
// tmp.put("bottomHeight",point.getBottomHeight());
// tmp.put("topHeight",point.getTopHeight());
// result.add(tmp);
// });
// });
//
// return result;
// }
/**
* 返回线要素的
* @param fillPointList
* @param radius
* @return
*/
private static SimpleFeature searchBourdary(List<Map<String, Object>> fillPointList, double radius,String name) {
AlphaShape alphaShape = new AlphaShape(radius);
// 构建顶点集合对象
List<Vector2D> points = new ArrayList<>();
fillPointList.forEach(pointMap->{
Point p = (Point)pointMap.get("geom");
Vector2D vector2D = new Vector2D(p.getX(), p.getY());
points.add(vector2D);
});
List<List<Vector2D>> lists = alphaShape.boundaryPoints(points);
List<Coordinate> coordList = new ArrayList();
lists.forEach(pointList->{
pointList.forEach(point->{
coordList.add(new Coordinate(point.getX(),point.getY()));
});
});
Coordinate[] CoordinateArray = coordList.stream().toArray(Coordinate[]::new);
LineString lineString = geometryFactory.createLineString(CoordinateArray);
// 构建添加线要素层
final SimpleFeatureType TYPE =createLineStringFeatureType();
SimpleFeatureBuilder featureBuilder = new SimpleFeatureBuilder(TYPE);
String uuid = UUID.randomUUID().toString();
featureBuilder.add(lineString);
featureBuilder.add(name);
featureBuilder.add(uuid);
SimpleFeature feature = featureBuilder.buildFeature(uuid);
// featureCollection.add(feature);
return feature;
}
private static Map<String, Object> calcGrid(List<List<Map<String, Object>>> pointRingFeaturesProperties, Point centroid,double meter) throws Exception {
List<Map<String,Object>> lastRing = pointRingFeaturesProperties.get(pointRingFeaturesProperties.size()-1);
// 定义地图单位的点距,创建缓冲区用的到
double length = isProjcetion?(meter/2):meter2degree(meter/2);
List<Map<String,Object>> fillPointList = new ArrayList<>();// 填方的点集
List<Map<String,Object>> cutPointList = new ArrayList<>();// 挖方的点集
double fillSum = 0;// 填方量
double cutSum = 0;// 挖方量
for(Map<String,Object> featuresProperties:lastRing){
Point point = (Point)featuresProperties.get("geom");
// 创建缓存面
LineString lineString = geometryFactory.createLineString(new Coordinate[]{new Coordinate(point.getX(),point.getY()),new Coordinate(centroid.getX(),centroid.getY())});
BufferParameters bufferParameters = new BufferParameters();
bufferParameters.setEndCapStyle(BufferParameters.CAP_FLAT);
bufferParameters.setJoinStyle(BufferParameters.JOIN_MITRE);
BufferOp bufOp = new BufferOp(lineString,bufferParameters);
Polygon bufferLinePolygon = (Polygon)bufOp.getResultGeometry(length);
// 计算
boolean isGoDown = true;// 是否继续往下
boolean isGoHead = true;
outer:
for(int j=0;j<pointRingFeaturesProperties.size()-1;j++){
List<Map<String,Object>> ringPointMap = pointRingFeaturesProperties.get(j);
for (Map<String, Object> pointMap : ringPointMap) {
Point p = (Point) pointMap.get("geom");
if (bufferLinePolygon.contains(p)) {
double bottomHeight = (double) pointMap.get("bottomHeight");
double topHeight = (double) pointMap.get("topHeight");
double currHeight = getElevation(p);
pointMap.put("elev",currHeight);
if (currHeight > bottomHeight && currHeight < topHeight) {
break outer;
}
if (currHeight >= topHeight && isGoHead) {
isGoDown = false;
if (!cutPointList.contains(pointMap)) {//去除重复点,避免重复计算
cutPointList.add(pointMap);
cutSum += (currHeight - topHeight) * Math.pow(meter, 2);
}
}
if (currHeight <= bottomHeight && isGoDown) {
isGoHead = false;
if (!fillPointList.contains(pointMap)) {//去除重复点,避免重复计算
fillPointList.add(pointMap);
fillSum += (currHeight - bottomHeight) * Math.pow(meter, 2);
}
}
}
}
}
}
Map<String, Object> result = new HashMap<>();
result.put("fillPointList",fillPointList);
result.put("cutPointList",cutPointList);
result.put("fillSum",fillSum);
result.put("cutSum",cutSum);
return result;
}
private static Map<String,Object> filterGrid(List<List<Map<String,Object>>> pointRingFeaturesProperties,GeometryFactory geometryFactory,Point centroid,double meter,String type){
Map<String,Object> result = new HashMap<>();
double sum = 0;
List<Map<String,Object>> filterGrid = new ArrayList<>();
for(int i=0;i<pointRingFeaturesProperties.size();i++){
List<Map<String,Object>> ringPointList = pointRingFeaturesProperties.get(i);
for(int j=0;j<ringPointList.size();j++){
Map<String,Object> pfp = ringPointList.get(j);
double height = (double)pfp.get("height");
Point p = (Point) pfp.get("geom");
try {
double elevation = getElevation(p);
// 填
if((type.equals("-1")&&height - elevation>=0)||(type.equals("1")&&height - elevation<0)){
List<Map<String, Object>> centerLinePointSet = getCenterLinePoint(geometryFactory,pointRingFeaturesProperties,i,centroid,meter,p);
if(centerLinePointSet.isEmpty()){
filterGrid.add(pfp);
sum += Math.abs(height - elevation)*Math.pow(meter, 2);
}else{
int count = 0;
for(int k=0;k<centerLinePointSet.size();k++){
Map<String, Object> lineOnPoint = centerLinePointSet.get(k);
double hp = (double)lineOnPoint.get("height");
Point geom = (Point) lineOnPoint.get("geom");
double ev = getElevation(geom);
if(ev>hp){
count++;
}
}
if(count==0){
filterGrid.add(pfp);
sum += Math.abs(height - elevation)*Math.pow(meter, 2);
}
}
}
} catch (Exception e) {
throw new RuntimeException(e);
}
}
};
result.put("sum",Integer.valueOf(type)*sum);
result.put("filterGrid",filterGrid);
return result;
}
private static List<Map<String,Object>> getCenterLinePoint( GeometryFactory geometryFactory,List<List<Map<String,Object>>> pointFeaturesProperties,int index,Point centroid,double side,Point point){
// 创建缓存面
LineString lineString = geometryFactory.createLineString(new Coordinate[]{new Coordinate(point.getX(),point.getY()),new Coordinate(centroid.getX(),centroid.getY())});
BufferParameters bufferParameters = new BufferParameters();
bufferParameters.setEndCapStyle(BufferParameters.CAP_FLAT);
bufferParameters.setJoinStyle(BufferParameters.JOIN_MITRE);
BufferOp bufOp = new BufferOp(lineString,bufferParameters);
double currDegree = side/2;
Polygon bufferLine = (Polygon)bufOp.getResultGeometry(currDegree);
List<Map<String, Object>> result = new ArrayList<>();
for(int i=0;i<index;i++){
List<Map<String,Object>> ring = pointFeaturesProperties.get(i);
for(int j=0;j<ring.size();j++){
Map<String, Object> pfp = ring.get(j);
Point p = (Point) pfp.get("geom");
if(bufferLine.intersects(p)){
result.add(pfp);
}
}
}
return result;
}
/**
* 创建矩形外扩的环点
* 打点方式 四条边顺时针按给定点距取点再buffer点距取矩形四条边重复之前按点距取点直到外扩的最后一个矩形
* @param rectangle 设计平面
* @param targetHeight 设计标高
* @param bottomRatio 往下放坡的坡度比
* @param topRatio 往上放坡的坡度比
* @param side 点距地图单位
* @return
*/
private static List<List<Map<String, Object>>> createRingPoint(Polygon rectangle, double targetHeight, double bottomRatio,double topRatio, double side) {
List<List<Map<String, Object>>> ringPoint = new ArrayList<>();
// 先在输入面上打点
List<Map<String, Object>> firstList = new ArrayList<>();
Coordinate[] coordinates = rectangle.getCoordinates();
Coordinate[] coordArr1 = new Coordinate[]{coordinates[0],coordinates[1]};
Coordinate[] coordArr2 = new Coordinate[]{coordinates[1],coordinates[2]};
Coordinate[] coordArr3 = new Coordinate[]{coordinates[2],coordinates[3]};
Coordinate[] coordArr4 = new Coordinate[]{coordinates[3],coordinates[0]};
firstList.addAll(addLinePoint(coordArr1,geometryFactory,side,targetHeight,targetHeight));
firstList.addAll(addLinePoint(coordArr2,geometryFactory,side,targetHeight,targetHeight));
firstList.addAll(addLinePoint(coordArr3,geometryFactory,side,targetHeight,targetHeight));
firstList.addAll(addLinePoint(coordArr4,geometryFactory,side,targetHeight,targetHeight));
ringPoint.add(firstList);
// 计算高程用的点距,单位为米
double length = isProjcetion?side:degree2meter(side);
for(int i=1;i<extendLen/length;i++){
List<Map<String, Object>> pointRingFeaturesProperties = new ArrayList<>();
double bottomHeight = targetHeight - length*bottomRatio*i;
double topHeight = targetHeight + length*topRatio*i;
BufferParameters bufferParameters = new BufferParameters();
bufferParameters.setEndCapStyle(BufferParameters.CAP_FLAT);
bufferParameters.setJoinStyle(BufferParameters.JOIN_MITRE);
BufferOp bufOp = new BufferOp(rectangle,bufferParameters);
double currDegree = side*i;
Geometry bg = bufOp.getResultGeometry(currDegree);// 缓冲分析的参数是地图单位的值
Coordinate[] boundCoord = bg.getBoundary().getCoordinates();
Coordinate[] boundArr1 = new Coordinate[]{boundCoord[0],boundCoord[1]};
Coordinate[] boundArr2 = new Coordinate[]{boundCoord[1],boundCoord[2]};
Coordinate[] boundArr3 = new Coordinate[]{boundCoord[2],boundCoord[3]};
Coordinate[] boundArr4 = new Coordinate[]{boundCoord[3],boundCoord[0]};
pointRingFeaturesProperties.addAll(addLinePoint(boundArr1,geometryFactory,side,bottomHeight,topHeight));
pointRingFeaturesProperties.addAll(addLinePoint(boundArr2,geometryFactory,side,bottomHeight,topHeight));
pointRingFeaturesProperties.addAll(addLinePoint(boundArr3,geometryFactory,side,bottomHeight,topHeight));
pointRingFeaturesProperties.addAll(addLinePoint(boundArr4,geometryFactory,side,bottomHeight,topHeight));
ringPoint.add(pointRingFeaturesProperties);
}
return ringPoint;
}
// 添加一个面的点集
public static List<List<Map<String,Object>>> addPartPoint(Coordinate[] bgcoords,
GeometryFactory geometryFactory,
double side,
Polygon p,
List<Map<String,Object>> pointFeaturesProperties,
double outHeight,
String type){
List<List<Map<String,Object>>> gridList = new ArrayList();
// 计算外边界上的点
List<Map<String,Object>> headPart = new ArrayList<>();
double t1 = (bgcoords[1].y - bgcoords[0].y)/(bgcoords[1].x-bgcoords[0].x);
Coordinate startCoord = new Coordinate(bgcoords[0].x,bgcoords[0].y);
Map<String,Object> centerPointMap = new HashMap<>();
Point center = geometryFactory.createPoint(startCoord);
centerPointMap.put("geom",center);
centerPointMap.put("id",UUID.randomUUID());
centerPointMap.put("bottomHeight",outHeight);
centerPointMap.put("topHeight",outHeight);
if(!type.equals("0")){
pointFeaturesProperties.add(centerPointMap);
headPart.add(centerPointMap);
}
List<Map<String,Object>> lastPart = addPointSet(bgcoords,geometryFactory,side,p,pointFeaturesProperties,center,outHeight,type);
headPart.addAll(lastPart);
gridList.add(headPart);
int len = (int)Math.ceil((Math.sqrt(Math.pow(bgcoords[1].y - bgcoords[0].y,2)+Math.pow(bgcoords[1].x - bgcoords[0].x,2)))/side);
//确定符号
double xx = bgcoords[1].x-bgcoords[0].x;
double yy = bgcoords[1].y-bgcoords[0].y;
int xSign = xx>=0?1:-1;
int ySign = yy>=0?1:-1;
for(int i=1;i<len;i++){
double nextX = startCoord.x + xSign * i*side*Math.abs(Math.cos(Math.atan(t1)));
double nextY = startCoord.y + ySign * i*side*Math.abs(Math.sin(Math.atan(t1)));
List<Map<String,Object>> extHeadPart = new ArrayList<>();
Coordinate nextCoord = new Coordinate(nextX,nextY);
Map<String,Object> nextPointMap = new HashMap<>();
Point nextPoint = geometryFactory.createPoint(nextCoord);
nextPointMap.put("geom",nextPoint);
nextPointMap.put("id",UUID.randomUUID());
nextPointMap.put("bottomHeight",outHeight);
nextPointMap.put("topHeight",outHeight);
if(type.equals("0")){
if(p.intersects(nextPoint)){
pointFeaturesProperties.add(nextPointMap);
extHeadPart.add(nextPointMap);
}
}else{
pointFeaturesProperties.add(nextPointMap);
extHeadPart.add(nextPointMap);
}
// 添加垂直于边界线的点阵
List<Map<String,Object>> extLastPart = addPointSet(bgcoords,geometryFactory,side,p,pointFeaturesProperties,nextPoint,outHeight,type);
extHeadPart.addAll(extLastPart);
gridList.add(extHeadPart);
}
return gridList;
}
/**
* 添加一个线的点集
* @param bgcoords 线的坐标
* @param geometryFactory
* @param side 点距地图单位
* @param bottomHeight 线所处的底部放坡高程
* @param topHeight 线做吃的顶部放坡高程
* @return 点集map数组
*/
public static List<Map<String,Object>> addLinePoint(Coordinate[] bgcoords, GeometryFactory geometryFactory, double side, double bottomHeight,double topHeight){
// double degree = side;
// 计算外边界上的点
List<Map<String,Object>> lineList = new ArrayList<>();
Coordinate startCoord = new Coordinate(bgcoords[0].x,bgcoords[0].y);
Map<String,Object> startPointMap = new HashMap<>();
Point sp = geometryFactory.createPoint(startCoord);
startPointMap.put("geom",sp);
startPointMap.put("id",UUID.randomUUID().toString());
startPointMap.put("bottomHeight",bottomHeight);
startPointMap.put("topHeight",topHeight);
lineList.add(startPointMap);
int len = (int)Math.floor((Math.sqrt(Math.pow(bgcoords[1].y - bgcoords[0].y,2)+Math.pow(bgcoords[1].x - bgcoords[0].x,2)))/side);
double t1 = (bgcoords[1].y - bgcoords[0].y)/(bgcoords[1].x-bgcoords[0].x);
//确定符号
double xx = bgcoords[1].x-bgcoords[0].x;
double yy = bgcoords[1].y-bgcoords[0].y;
int xSign = xx>=0?1:-1;
int ySign = yy>=0?1:-1;
// 计算坐标直接用地图单位的点距
for(int i=1;i<len;i++){
double nextX = startCoord.x + xSign * i*side*Math.abs(Math.cos(Math.atan(t1)));
double nextY = startCoord.y + ySign * i*side*Math.abs(Math.sin(Math.atan(t1)));
if(xx == 0){
nextX = startCoord.x;
nextY = startCoord.y + ySign * i*side;
}
Coordinate nextCoord = new Coordinate(nextX,nextY);
Map<String,Object> nextPointMap = new HashMap<>();
Point nextPoint = geometryFactory.createPoint(nextCoord);
nextPointMap.put("geom",nextPoint);
nextPointMap.put("id",UUID.randomUUID().toString());
nextPointMap.put("bottomHeight",bottomHeight);
nextPointMap.put("topHeight",topHeight);
lineList.add(nextPointMap);
}
return lineList;
}
public static List<Map<String,Object>> addPointSet(Coordinate[] bgcoords,
GeometryFactory geometryFactory,
double side,
Polygon polygon,
List<Map<String,Object>> pointFeaturesProperties,
Point startPoint,
double boundHeight,
String type){
List<Map<String,Object>> extendPointProperties = new ArrayList<>();
double t1 = (bgcoords[1].y - bgcoords[0].y)/(bgcoords[1].x-bgcoords[0].x);
//确定符号
double xx = bgcoords[1].x-bgcoords[0].x;
double yy = bgcoords[1].y-bgcoords[0].y;
int xSign = (xx>=0&&yy>0)||(xx<0&&yy>0)?1:-1;
int ySign = (xx<0&&yy<0)||(xx<0&&yy>0)?1:-1;
// 算高程要用到,地图单位到米
double length = isProjcetion?side:degree2meter(side);
for(int i=1;i<3000;i++){
double x1 = xSign * Math.abs(Math.sin(Math.atan(t1)))*side*i + startPoint.getX();
double y1 = ySign * Math.abs(Math.cos(Math.atan(t1)))*side*i + startPoint.getY();
// 当垂直与坐标轴时
if(bgcoords[1].x==bgcoords[0].x){
x1 = side*i + startPoint.getX();
y1 = startPoint.getY();
}
Point point = geometryFactory.createPoint(new Coordinate(x1,y1));
boolean contains = polygon.contains(point);
if(contains){
Map<String,Object> tmpPoint = new HashMap<>();
tmpPoint.put("geom",point);
tmpPoint.put("id",UUID.randomUUID());
double ph = boundHeight;//+(degree2meter(side))*i/2;
if(type.equals("-1")){//下挖的高程计算
ph+=(length)*i/2;
}else if(type.equals("1")){//
ph-=(length)*i/1.5;
}
tmpPoint.put("bottomHeight",ph);
tmpPoint.put("topHeight",ph);
pointFeaturesProperties.add(tmpPoint);
extendPointProperties.add(tmpPoint);
}
}
return extendPointProperties;
}
public static void loadDem(String path) throws Exception{
File file = new File(path);
// File file = new File("src\\main\\resources\\public\\data\\CSDEM.tif");
// File file = new File("E:\\data\\CB测试\\sample_N30E119.tif");
// File file = new File("E:\\data\\CB测试\\sample_N30E119.tif");
// Hints.FORCE_LONGITUDE_FIRST_AXIS_ORDER设置经度为第一轴顺序
GeoTiffReader reader = new GeoTiffReader(file, new Hints(Hints.FORCE_LONGITUDE_FIRST_AXIS_ORDER, Boolean.TRUE));
// GridCoverage2D coverage = reader.read(null);
// //设置坐标系
// CoordinateReferenceSystem targetCRS = CRS.decode("EPSG:4326");
// // 将 GridCoverage 进行重采样转换为另一个 CRS
// CreateRingGridFixProjection.coverage = (GridCoverage2D) Operations.DEFAULT.resample(coverage, targetCRS);
CreateRingGridFixGeo.coverage = reader.read(null);
}
public static double getElevation(Point p) throws Exception{
// 设置经纬度及坐标系等信息
Position position = new Position2D(p.getX(),p.getY());
GridGeometry2D gridGeometry = CreateRingGridFixGeo.coverage.getGridGeometry();
GridCoordinates2D gridPoint = gridGeometry.worldToGrid(position);
double[] elevation = new double[1];
CreateRingGridFixGeo.coverage.evaluate(gridPoint,elevation);
double height = elevation[0]; // 这里的高度就是查询位置的高程值
// System.out.println("---------"+height);//---------108.0
return height;
}
public static double meter2degree(double meter){
return meter/(Math.PI*6371004/180);
}
public static double degree2meter(double degree){
return degree*(Math.PI*6371004)/180;
}
/**
* 对栅格重采样
* @param envelope 限定的矩形几何
* @param inputPath 输入的栅格路径
* @param outputPath 输出的栅格路径
* @param r1 东西方向像元大小地图单位
* @param r2 南北方向像元大小地图单位
*/
public static void reSample(Geometry envelope,String inputPath, String outputPath, double r1, double r2) {
// 打开输入栅格文件
Dataset dataset = gdal.Open(inputPath, gdalconst.GA_ReadOnly);
double[] gt1 = dataset.GetGeoTransform();
Coordinate[] envCoords = envelope.getCoordinates();
double envWidth = (envCoords[2].x - envCoords[0].x);
double envHeight = (envCoords[2].y - envCoords[0].y);
int xSize = new Double(envWidth/ r1).intValue();
int ySize = new Double(envHeight / r2).intValue();
ySize = Math.abs(ySize);
double[] gt2 = {envCoords[0].x, r1, 0.0, envCoords[0].y, 0.0, r2};
Driver driver = gdal.GetDriverByName("GTiff");
Dataset datasetWarp = driver.Create(outputPath, xSize, ySize, gdalconst.GA_Update, gdalconst.GDT_Float32);
datasetWarp.SetGeoTransform(gt2);
datasetWarp.SetSpatialRef(dataset.GetSpatialRef());
Vector<String> vector = new Vector<>();
vector.add("-r");
vector.add("cubicspline");//执行三次卷积插值法重采样
WarpOptions warpOptions = new WarpOptions(vector);
gdal.Warp(datasetWarp, new Dataset[]{dataset}, warpOptions);
datasetWarp.FlushCache();
datasetWarp.delete();
dataset.delete();
}
}