diff --git a/pom.xml b/pom.xml
index a75bc57..914536c 100644
--- a/pom.xml
+++ b/pom.xml
@@ -93,6 +93,11 @@
gt-opengis
29.5
+
+ org.geotools
+ gt-geotiff
+ ${geotools.version}
+
org.gdal
diff --git a/src/main/java/org/example/test/CreateContour.java b/src/main/java/org/example/test/CreateContour.java
new file mode 100644
index 0000000..13bda75
--- /dev/null
+++ b/src/main/java/org/example/test/CreateContour.java
@@ -0,0 +1,288 @@
+package org.example.test;
+
+import it.geosolutions.jaiext.range.Range;
+import org.gdal.gdal.Band;
+import org.gdal.gdal.Dataset;
+import org.gdal.gdal.gdal;
+import org.gdal.gdalconst.gdalconst;
+import org.gdal.ogr.Layer;
+import org.gdal.osr.SpatialReference;
+import org.geotools.api.data.SimpleFeatureSource;
+import org.geotools.api.data.SimpleFeatureStore;
+import org.geotools.api.data.Transaction;
+import org.geotools.api.feature.simple.SimpleFeature;
+import org.geotools.api.feature.simple.SimpleFeatureType;
+import org.geotools.api.geometry.Position;
+import org.geotools.api.referencing.crs.CoordinateReferenceSystem;
+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.DefaultTransaction;
+import org.geotools.data.collection.ListFeatureCollection;
+import org.geotools.data.shapefile.ShapefileDataStore;
+import org.geotools.data.simple.SimpleFeatureCollection;
+import org.geotools.feature.simple.SimpleFeatureBuilder;
+import org.geotools.feature.simple.SimpleFeatureTypeBuilder;
+import org.geotools.gce.geotiff.GeoTiffReader;
+import org.geotools.geometry.Position2D;
+import org.geotools.geometry.jts.ReferencedEnvelope;
+import org.geotools.referencing.CRS;
+import org.geotools.referencing.crs.DefaultGeographicCRS;
+import org.geotools.util.factory.Hints;
+import org.locationtech.jts.geom.*;
+import org.geotools.data.shapefile.ShapefileDataStoreFactory;
+import org.locationtech.jts.operation.buffer.BufferOp;
+import org.locationtech.jts.operation.buffer.BufferParameters;
+import org.opengis.geometry.Envelope;
+
+import java.io.File;
+import java.io.IOException;
+import java.io.Serializable;
+import java.nio.charset.Charset;
+import java.util.*;
+
+public class CreateContour {
+ public static void main(String[] args) throws Exception{
+ loadDem();
+ // 假设我们使用WGS84坐标系统
+ CoordinateReferenceSystem crs = CRS.decode("EPSG:4326");
+ GeometryFactory geometryFactory = new GeometryFactory();
+ // 定义矩形的左下角和右上角坐标(WGS84)
+ List coords = new ArrayList(){{
+ add(new double[]{113.52896835361318, 23.654018962759377});
+ add(new double[]{113.529270627066524, 23.654027766840539 });
+ add(new double[]{ 113.529294104616298, 23.653772448486752});
+ add(new double[]{ 113.528977157694356, 23.653763644405586 });
+ add(new double[]{113.52896835361318, 23.654018962759377});
+ }};
+ // 创建矩形的四个角点坐标
+ 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;
+
+ System.out.println("设计标高:"+targetHeight);
+
+ double allSum = calcBlanceHeight(coordinates,geometryFactory,targetHeight);
+// double[] stepSet = new double[]{5,3,1,0.5,0.25,0.125,0.0625,0.03125,0.0156125,0.007806125};
+// PrimitiveIterator.OfDouble iterator = Arrays.stream(stepSet).iterator();
+// double nextHeight=targetHeight,step = iterator.next();
+// while (Math.abs(allSum)-100>0){
+// nextHeight = allSum>=0?nextHeight-step:nextHeight+step;
+// System.out.println("nextHeight:"+nextHeight);
+// double firstSum = calcBlanceHeight(coordinates,geometryFactory,nextHeight);
+// if(Math.abs(firstSum)-1<0){
+// break;
+// }
+// nextHeight = firstSum>=0?nextHeight-step:nextHeight+step;
+// System.out.println("nextHeight:"+nextHeight);
+// System.out.println("firstSum:"+firstSum);
+// allSum = calcBlanceHeight(coordinates,geometryFactory,nextHeight);
+// if((firstSum>0&&allSum<0)||(firstSum<0&&allSum>0)){
+// step = iterator.next();
+// }
+// System.out.println("allSum:"+allSum);
+// }
+// System.out.println("lastHeight:"+nextHeight);
+
+ }
+ public static double calcBlanceHeight(Coordinate[] coordinates,GeometryFactory geometryFactory,double height) throws Exception{
+ double bufferLen = 100, stepLen = 2;
+
+ List featuresList = new ArrayList<>();
+ SimpleFeatureType TYPE = createFeatureType();
+ SimpleFeatureBuilder featureBuilder = new SimpleFeatureBuilder(TYPE);
+
+ Polygon rectangle = geometryFactory.createPolygon(coordinates);
+ LineString boundary = (LineString)rectangle.getBoundary();
+ featureBuilder.add(geometryFactory.createLineString(boundary.getCoordinates()));
+ featureBuilder.add("目标平面的第一条等高线");
+ featureBuilder.add(height);
+ //buildFeature(id)创建一个要素,ID可以是null,当ID为null时会由builder内部生成
+ SimpleFeature feature = featureBuilder.buildFeature("1");
+ featuresList.add(feature);
+
+ // 计算放坡的等高线
+ for(int i=1;i<=bufferLen/stepLen;i++){
+ BufferParameters bufferParameters = new BufferParameters();
+ bufferParameters.setEndCapStyle(BufferParameters.CAP_FLAT);
+ bufferParameters.setJoinStyle(BufferParameters.JOIN_MITRE);
+ BufferOp bufOp = new BufferOp(rectangle,bufferParameters);
+
+ double degree = meter2degree(stepLen*i);
+ Polygon bg = (Polygon)bufOp.getResultGeometry(degree);
+ LineString lineString = (LineString)bg.getBoundary();
+ featureBuilder.add(lineString);
+ featureBuilder.add(("扩展平面的第"+i+"条等高线"));
+ featureBuilder.add(height-i*stepLen/2);
+ //buildFeature(id)创建一个要素,ID可以是null,当ID为null时会由builder内部生成
+ SimpleFeature simpleFeature = featureBuilder.buildFeature("ext"+i);
+ featuresList.add(simpleFeature);
+ }
+ SimpleFeatureCollection collection = new ListFeatureCollection(TYPE, featuresList);
+ ReferencedEnvelope bounds = collection.getBounds();
+
+ ShapefileDataStore newDataStore = getShapefileDataStore(collection);
+ newDataStore.createSchema(TYPE);
+ newDataStore.setCharset(Charset.forName("UTF-8"));
+ Transaction transaction = new DefaultTransaction("create");
+
+ String typeName = newDataStore.getTypeNames()[0];
+ SimpleFeatureSource featureSource = newDataStore.getFeatureSource(typeName);
+ SimpleFeatureType SHAPE_TYPE = featureSource.getSchema();
+ /*
+ * The Shapefile format has a couple limitations:
+ * - "the_geom" is always first, and used for the geometry attribute name
+ * - "the_geom" must be of type Point, MultiPoint, MuiltiLineString, MultiPolygon
+ * - Attribute names are limited in length
+ * - Not all data types are supported (example Timestamp represented as Date)
+ *
+ * Each data store has different limitations so check the resulting SimpleFeatureType.
+ */
+ System.out.println("SHAPE:" + SHAPE_TYPE);
+
+ if (featureSource instanceof SimpleFeatureStore) {
+ SimpleFeatureStore featureStore = (SimpleFeatureStore) featureSource;
+ /*
+ * SimpleFeatureStore has a method to add features from a
+ * SimpleFeatureCollection object, so we use the ListFeatureCollection
+ * class to wrap our list of features.
+ */
+
+ featureStore.setTransaction(transaction);
+ try {
+ featureStore.addFeatures(collection);
+ transaction.commit();
+ } catch (Exception problem) {
+ problem.printStackTrace();
+ transaction.rollback();
+ } finally {
+ transaction.close();
+ }
+// System.exit(0); // success!
+ } else {
+ System.out.println(typeName + " does not support read/write access");
+ System.exit(1);
+ }
+
+ // 注册所有gdal的驱动
+ gdal.AllRegister();
+ Dataset vector = gdal.OpenEx("E:/code/CutFillDemo/target/classes/tmp.shp",gdalconst.OF_VECTOR,null,null,null);
+ int layerCount = vector.GetLayerCount();
+ if(layerCount == 0){
+ System.err.println("矢量数据没有图层!");
+ System.exit(1);
+ }
+ Layer layer = vector.GetLayerByIndex(0);
+ SpatialReference crs = layer.GetSpatialRef();
+ double[] boundsArr = layer.GetExtent();
+ double pixelSize = 0.00025;
+ int NoData_value = -9999;
+ int x_res = (int)((boundsArr[1] - boundsArr[0]) / pixelSize);
+ int y_res = (int)((boundsArr[3] - boundsArr[2]) / pixelSize);
+ if(x_res == 0){
+ x_res = 2000;
+ }
+ if(y_res == 0){
+ y_res = 2000;
+ }
+ // 创建dem输出的数据集
+ String demOutPath = "E:\\data\\demOutFile.tiff";
+ Dataset target_ds = gdal.GetDriverByName("GTiff").Create(demOutPath,x_res,y_res,1, gdalconst.GDT_Byte);
+ target_ds.SetGeoTransform(new double[]{ boundsArr[0], pixelSize,0,boundsArr[3],0,-pixelSize });
+ target_ds.SetProjection(crs.ExportToWkt());
+
+
+ Band band = target_ds.GetRasterBand(1);
+ band.SetNoDataValue(NoData_value);
+// band.FlushCache();
+ Vector