完成打点,添加填挖计算

This commit is contained in:
user 2024-04-18 17:32:30 +08:00
parent 1307698f4a
commit 9bef4281a4
7 changed files with 224 additions and 147 deletions

10
pom.xml
View File

@ -82,5 +82,15 @@
<artifactId>gt-grid</artifactId>
<version>${geotools.version}</version>
</dependency>
<dependency>
<groupId>org.geotools</groupId>
<artifactId>gt-process</artifactId>
<version>${geotools.version}</version>
</dependency>
<dependency>
<groupId>org.geotools</groupId>
<artifactId>gt-opengis</artifactId>
<version>29.5</version>
</dependency>
</dependencies>
</project>

View File

@ -26,6 +26,7 @@ public class Caclelate {
* 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);

View File

@ -1,9 +1,6 @@
package org.example.test;
import org.example.utils.IntersectionBuilder;
import org.geotools.api.data.FeatureSource;
import org.geotools.api.data.FileDataStore;
import org.geotools.api.data.FileDataStoreFinder;
import org.geotools.api.data.SimpleFeatureSource;
import org.geotools.api.feature.simple.SimpleFeature;
import org.geotools.api.feature.simple.SimpleFeatureType;
@ -14,39 +11,29 @@ 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.data.DataUtilities;
import org.geotools.data.collection.CollectionFeatureSource;
import org.geotools.data.collection.ListFeatureCollection;
import org.geotools.data.simple.SimpleFeatureCollection;
import org.geotools.data.store.EmptyFeatureCollection;
import org.geotools.feature.DefaultFeatureCollection;
import org.geotools.feature.FeatureCollection;
import org.geotools.feature.simple.SimpleFeatureBuilder;
import org.geotools.feature.simple.SimpleFeatureTypeBuilder;
import org.geotools.filter.FilterFactoryImpl;
import org.geotools.filter.function.color.ConstrastFunction;
import org.geotools.gce.geotiff.GeoTiffReader;
import org.geotools.geometry.Position2D;
import org.geotools.geometry.jts.JTSFactoryFinder;
import org.geotools.geometry.jts.ReferencedEnvelope;
import org.geotools.grid.Grids;
import org.geotools.map.FeatureLayer;
import org.geotools.map.Layer;
import org.geotools.map.MapContent;
import org.geotools.referencing.CRS;
import org.geotools.referencing.crs.DefaultGeographicCRS;
import org.geotools.styling.FillImpl;
import org.geotools.styling.SLD;
import org.geotools.swing.JMapFrame;
import org.geotools.swing.data.JFileDataStoreChooser;
import org.geotools.util.factory.Hints;
import org.locationtech.jts.geom.*;
import org.locationtech.jts.operation.buffer.BufferOp;
import org.locationtech.jts.operation.buffer.BufferParameters;
import java.awt.image.RenderedImage;
import java.io.File;
import java.util.*;
import java.util.concurrent.atomic.AtomicReference;
public class CreateGrid {
public static GridCoverage2D coverage;
@ -76,6 +63,7 @@ public class CreateGrid {
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("height", double.class);
return builder.buildFeatureType();
}
public static void main(String[] args) throws Exception{
@ -99,16 +87,16 @@ public class CreateGrid {
coordinates[3] = new Coordinate(coords.get(3)[0], coords.get(3)[1]);
// 闭合多边形复制第一个坐标作为最后一个坐标
coordinates[4] = new Coordinate(coords.get(0)[0], coords.get(0)[1]);
// coordinates[0] = new Coordinate(0, 0);
// coordinates[1] = new Coordinate(0, 30);
// coordinates[2] = new Coordinate(30, 30);
// coordinates[3] = new Coordinate(30, 0);
// // 闭合多边形复制第一个坐标作为最后一个坐标
// coordinates[4] = new Coordinate(0, 0);
// 创建GeometryFactory
GeometryFactory geometryFactory = new GeometryFactory();
List<Map<String,Object>> featuresProperties = new ArrayList<>();
List<Map<String,Object>> pointFeaturesProperties = new ArrayList<>();
List<Map<String,Object>> featuresProperties = new ArrayList<>();// 面要素属性集合
// List<Map<String,Object>> pointFeaturesProperties = new ArrayList<>();// 扩展点属性集合
List<Map<String,Object>> tPointFeaturesProperties = new ArrayList<>();// 目标平面的点集属性集合
List<Map<String,Object>> p1PointFeaturesProperties = new ArrayList<>();// p1面内的扩展点属性集合
List<Map<String,Object>> p2PointFeaturesProperties = new ArrayList<>();// p2面内的扩展点属性集合
List<Map<String,Object>> p3PointFeaturesProperties = new ArrayList<>();// p3面内的扩展点属性集合
List<Map<String,Object>> p4PointFeaturesProperties = new ArrayList<>();// p4面内的扩展点属性集合
// 创建多边形矩形
Map<String,Object> innerPolygon = new HashMap<>();
@ -118,13 +106,20 @@ public class CreateGrid {
innerPolygon.put("id",1);
featuresProperties.add(innerPolygon);
// 获取中心点的高程
double targetHeight = getElevation(rectangle.getCentroid());
// 计算设计标高取四个角的高程平均值
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;
System.out.println("设计标高:"+targetHeight);
// 缓冲区建立外扩100m
double meter = 100;
double degree = meter/(2*Math.PI*6371004)*360;
double degree = meter2degree(meter);
double side = meter2degree(2);// side为间隔单位为米
BufferParameters bufferParameters = new BufferParameters();
bufferParameters.setEndCapStyle(BufferParameters.CAP_FLAT);
@ -150,7 +145,7 @@ public class CreateGrid {
extPolygon2.put("geom",p2);
extPolygon2.put("name","第二个面");
extPolygon2.put("id",3);
// featuresProperties.add(extPolygon2);
featuresProperties.add(extPolygon2);
//
Map<String,Object> extPolygon3 = new HashMap<>();
Coordinate[] cords3 = new Coordinate[] {bgcoords[2],bgcoords[3],pcoords[3],pcoords[2],bgcoords[2]};
@ -158,7 +153,7 @@ public class CreateGrid {
extPolygon3.put("geom",p3);
extPolygon3.put("name","第3个面");
extPolygon3.put("id",4);
// featuresProperties.add(extPolygon3);
featuresProperties.add(extPolygon3);
Map<String,Object> extPolygon4 = new HashMap<>();
Coordinate[] cords4 = new Coordinate[] {bgcoords[3],bgcoords[0],pcoords[0],pcoords[3],bgcoords[3]};
@ -166,20 +161,43 @@ public class CreateGrid {
extPolygon4.put("geom",p4);
extPolygon4.put("name","第4个面");
extPolygon4.put("id",5);
// featuresProperties.add(extPolygon4);
featuresProperties.add(extPolygon4);
// 外接矩形打点
Map<String,Object> outPolygon1 = new HashMap<>();
// Coordinate[] cords4 = new Coordinate[] {bgcoords[3],bgcoords[0],pcoords[0],pcoords[3],bgcoords[3]};
Polygon out1 = (Polygon) p1.getEnvelope();//geometryFactory.createPolygon(cords4);
Polygon out1 = (Polygon) rectangle.getEnvelope();//geometryFactory.createPolygon(cords4);
outPolygon1.put("geom",out1);
outPolygon1.put("name","第1个面的外边框");
outPolygon1.put("name","输入面的外接矩形");
outPolygon1.put("id",6);
// featuresProperties.add(outPolygon1);
featuresProperties.add(outPolygon1);
Coordinate[] outCoords = out1.getCoordinates();
Coordinate[] outLine = new Coordinate[]{outCoords[0],outCoords[1]};
addPartPoint(outLine,geometryFactory,side,rectangle,tPointFeaturesProperties,targetHeight,"0");
// 打点
double outBottomHeight = targetHeight - meter/2;
System.out.println("边界高程:"+outBottomHeight);
Coordinate[] coordArr1 = new Coordinate[]{bgcoords[0],bgcoords[1]};
List gridsP1 = addPartPoint(coordArr1,geometryFactory,side,p1,p1PointFeaturesProperties,outBottomHeight,"-1");
Coordinate[] coordArr2 = new Coordinate[]{bgcoords[1],bgcoords[2]};
List gridsP2 = addPartPoint(coordArr2,geometryFactory,side,p2,p2PointFeaturesProperties,outBottomHeight,"-1");
Coordinate[] coordArr3 = new Coordinate[]{bgcoords[2],bgcoords[3]};
List gridsP3 = addPartPoint(coordArr3,geometryFactory,side,p3,p3PointFeaturesProperties,outBottomHeight,"-1");
Coordinate[] coordArr4 = new Coordinate[]{bgcoords[3],bgcoords[0]};
List gridsP4 = addPartPoint(coordArr4,geometryFactory,side,p4,p4PointFeaturesProperties,outBottomHeight,"-1");
extPanelCalc(gridsP2,side);
// 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"));
@ -187,12 +205,39 @@ public class CreateGrid {
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(pointFeaturesProperties);
add(tPointFeaturesProperties);
add(p1PointFeaturesProperties);
add(p2PointFeaturesProperties);
add(p3PointFeaturesProperties);
add(p4PointFeaturesProperties);
}};
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("height"));
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);
}
});
// ReferencedEnvelope bounds = featureCollection.getBounds();
// 创建grids点side为间隔单位为米
double side = 1/(2*Math.PI*6371004)*360;
// 添加网格图层
// 创建grids点
ListFeatureCollection collection = new ListFeatureCollection(TYPE);
Iterator<SimpleFeature> iterator = featureCollection.stream().iterator();
iterator.next();
@ -200,137 +245,151 @@ public class CreateGrid {
SimpleFeatureSource source = new CollectionFeatureSource(collection);
IntersectionBuilder intersectionBuilder = new IntersectionBuilder(TYPE, source);
SimpleFeatureSource grid = Grids.createSquareGrid(collection.getBounds(), side, -1);
// 打点
Coordinate pc1 = bgcoords[0];
Coordinate pc2 = bgcoords[1];
Coordinate pc3 = bgcoords[2];
double t1 = (pc2.y - pc1.y)/(pc2.x-pc1.x);
double t2 = (pc3.y - pc2.y)/(pc3.x-pc2.x);
System.out.println("t1:"+t1);
System.out.println("t2:"+t2);
List<Geometry> centerPointList = new ArrayList();
Coordinate center1 = new Coordinate((pc1.x+pc2.x)/2,(pc1.y+pc2.y)/2);
Map<String,Object> centerPointMap = new HashMap<>();
Point center = geometryFactory.createPoint(center1);
// centerPointList.add(center);
// centerPointMap.put("geom",center);
// centerPointMap.put("id",10001);
// pointFeaturesProperties.add(centerPointMap);
//
// for(int i=1;i<200;i++){
// double x1 = Math.sin(Math.atan(t1))*side*i + center1.x;
// double y1 = - Math.cos(Math.atan(t1))*side*i + center1.y;
//
// Point p = geometryFactory.createPoint(new Coordinate(x1,y1));
//// double height = getElevation(p);
// boolean contains = p1.contains(p);
// if(contains){
// Map<String,Object> tmpPoint = new HashMap<>();
// centerPointList.add(p);
// tmpPoint.put("geom",p);
// tmpPoint.put("id",i);
// pointFeaturesProperties.add(tmpPoint);
// }else{
// break;
// }
// }
// 下一个点
// double nextX= center1.x + side*Math.cos(Math.atan(t1));
// double nextY = center1.y + side*Math.sin(Math.atan(t1));
// for(int i=1;i<200;i++){
// double x1 = Math.sin(Math.atan(t1))*side*i + nextX;
// double y1 = - Math.cos(Math.atan(t1))*side*i + nextY;
//
// Point p = geometryFactory.createPoint(new Coordinate(x1,y1));
//// double height = getElevation(p);
// boolean contains = p1.contains(p);
// if(contains){
// Map<String,Object> tmpPoint = new HashMap<>();
// centerPointList.add(p);
// tmpPoint.put("geom",p);
// tmpPoint.put("id",2000+i);
// pointFeaturesProperties.add(tmpPoint);
// }else{
// break;
// }
// }
Coordinate[] coordArr = new Coordinate[]{bgcoords[0],bgcoords[1],bgcoords[2]};
addPointSet(coordArr,geometryFactory,side,p1,pointFeaturesProperties,center);
final SimpleFeatureType PTYPE =createPointFeatureType();
DefaultFeatureCollection pointFeatureCollection = new DefaultFeatureCollection();
SimpleFeatureBuilder pointFeatureBuilder = new SimpleFeatureBuilder(PTYPE);
pointFeaturesProperties.forEach(props->{
pointFeatureBuilder.add(props.get("geom"));
pointFeatureBuilder.add(props.get("id"));
SimpleFeature feature = pointFeatureBuilder.buildFeature((String.valueOf(props.get("id"))));
pointFeatureCollection.add(feature);
});
// Create a map content and add our shapefile to it
MapContent map = new MapContent();
map.setTitle("填挖方计算-创建网格");
// 添加初始范围和扩展面图层
Style style = SLD.createSimpleStyle(featureCollection.getSchema());
Layer layer = new FeatureLayer(featureCollection, style);
map.addLayer(layer);
// 添加点
Style pointStyle = SLD.createSimpleStyle(pointFeatureCollection.getSchema());
Layer pointLlayer = new FeatureLayer(pointFeatureCollection, pointStyle);
map.addLayer(pointLlayer);
// 添加网格图层
Style styleBox = SLD.createSimpleStyle(grid.getSchema());
Layer layerBox = new FeatureLayer(grid, styleBox);
// map.addLayer(layerBox);
// Now display the map
JMapFrame.showMap(map);
targetPanelCalc(tPointFeaturesProperties,side);
}
public static void addPointSet(Coordinate[] bgcoords,
public static void targetPanelCalc (List<Map<String,Object>> pointFeaturesProperties,double side){
AtomicReference<Double> sum = new AtomicReference<>((double) 0);// 土方量总和负数挖整数填
pointFeaturesProperties.stream().forEach(pfp->{
double h = (double)pfp.get("height");
Point p = (Point)pfp.get("geom");
try {
double elevation = getElevation(p);
sum.updateAndGet(v -> new Double((double) (v + (h - elevation) * Math.pow(degree2meter(side), 2))));
} catch (Exception e) {
throw new RuntimeException(e);
}
});
System.out.println("sum:"+sum);
}
public static void extPanelCalc(List<List<Map<String,Object>>> gridPolygon,double side){
AtomicReference<Double> sum = new AtomicReference<>((double) 0);// 土方量总和整数填
for(int i=0;i<gridPolygon.size();i++){
List<Map<String,Object>> col = gridPolygon.get(i);
// 逆序比较
for(int j=col.size()-1;j>=0;j--){
Map<String,Object> pfp = col.get(j);
double h = (double)pfp.get("height");
Point p = (Point)pfp.get("geom");
try {
double elevation = getElevation(p);
if(elevation>h){
break;
}
sum.updateAndGet(v -> new Double((double) (v + (h - elevation) * Math.pow(degree2meter(side), 2))));
} catch (Exception e) {
throw new RuntimeException(e);
}
}
}
System.out.println("gridPolygon sum:"+sum);
}
// 添加一个面的点集
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("height",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("height",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;
}
public static List<Map<String,Object>> addPointSet(Coordinate[] bgcoords,
GeometryFactory geometryFactory,
double side,
Polygon p1,
Polygon polygon,
List<Map<String,Object>> pointFeaturesProperties,
Point startPoint){
Coordinate pc1 = bgcoords[0];
Coordinate pc2 = bgcoords[1];
Coordinate pc3 = bgcoords[2];
Point startPoint,
double boundHeight,
String type){
List<Map<String,Object>> extendPointProperties = new ArrayList<>();
double t1 = (pc2.y - pc1.y)/(pc2.x-pc1.x);
double t2 = (pc3.y - pc2.y)/(pc3.x-pc2.x);
System.out.println("t1:"+t1);
System.out.println("t2:"+t2);
List<Geometry> centerPointList = new ArrayList();
Coordinate center1 = new Coordinate((pc1.x+pc2.x)/2,(pc1.y+pc2.y)/2);
Map<String,Object> centerPointMap = new HashMap<>();
Point center = geometryFactory.createPoint(center1);
centerPointList.add(center);
centerPointMap.put("geom",startPoint);
centerPointMap.put("id",UUID.randomUUID());
pointFeaturesProperties.add(centerPointMap);
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;
for(int i=1;i<200;i++){
double x1 = Math.sin(Math.atan(t1))*side*i + startPoint.getX();
double y1 = - Math.cos(Math.atan(t1))*side*i + startPoint.getY();
Point p = geometryFactory.createPoint(new Coordinate(x1,y1));
boolean contains = p1.contains(p);
for(int i=1;i<300;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();
Point point = geometryFactory.createPoint(new Coordinate(x1,y1));
boolean contains = polygon.contains(point);
if(contains){
Map<String,Object> tmpPoint = new HashMap<>();
centerPointList.add(p);
tmpPoint.put("geom",p);
tmpPoint.put("geom",point);
tmpPoint.put("id",UUID.randomUUID());
double ph = boundHeight;//+(degree2meter(side))*i/2;
if(type.equals("-1")){//下挖的高程计算
ph+=(degree2meter(side))*i/2;
}else if(type.equals("1")){//
// ph-=(degree2meter(side))*i/2;
}
tmpPoint.put("height",ph);
pointFeaturesProperties.add(tmpPoint);
extendPointProperties.add(tmpPoint);
}else{
break;
}
}
return extendPointProperties;
}
public static void loadDem() throws Exception{
File file = new File("src\\main\\resources\\public\\data\\ASTGTMV003_N23E113_dem.tif");
@ -353,4 +412,11 @@ public class CreateGrid {
// 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;
}
}

Binary file not shown.