feat: 实现手势调整坡度

app/build.gradle

@@ -1,6 +1,7 @@
 plugins {
     alias(libs.plugins.android.application)
     alias(libs.plugins.kotlin.android)
+    alias(libs.plugins.kotlin.serialization)
 }
 
 android {
@@ -23,12 +24,15 @@ android {
             proguardFiles getDefaultProguardFile('proguard-android-optimize.txt'), 'proguard-rules.pro'
         }
     }
+    buildFeatures {
+        viewBinding true
+    }
     compileOptions {
-        sourceCompatibility JavaVersion.VERSION_11
-        targetCompatibility JavaVersion.VERSION_11
+        sourceCompatibility JavaVersion.VERSION_17
+        targetCompatibility JavaVersion.VERSION_17
     }
     kotlinOptions {
-        jvmTarget = '11'
+        jvmTarget = '17'
     }
 }
 
@@ -42,6 +46,15 @@ dependencies {
     implementation libs.androidx.constraintlayout
     implementation project(':delaunator')
     implementation project(':math')
+    implementation 'org.gdal:gdal:3.12.0'
+    implementation libs.kotlinx.serialization.json
+    implementation libs.ktor.client.core
+    implementation libs.ktor.client.cio
+    implementation libs.ktor.serialization.kotlinx.json
+    implementation libs.ktor.client.content.negotiation
+    implementation libs.ktor.client.logging
+    // https://mvnrepository.com/artifact/org.openrndr.extra/orx-marching-squares-jvm
+    implementation 'org.openrndr.extra:orx-marching-squares-jvm:0.4.5-alpha6'
 
     testImplementation libs.junit
     androidTestImplementation libs.androidx.junit

app/src/main/java/com/icegps/geotools/ContourConnector.kt

@@ -0,0 +1,402 @@
+package com.icegps.geotools
+
+import android.util.Log
+import com.icegps.common.helper.GeoHelper
+import com.icegps.geotools.ktx.TAG
+import com.icegps.geotools.ktx.niceStr
+import com.icegps.math.geometry.Vector3D
+import com.mapbox.geojson.Feature
+import com.mapbox.geojson.LineString
+import com.mapbox.geojson.Point
+import com.mapbox.maps.MapView
+import kotlin.math.abs
+import kotlin.math.floor
+import kotlin.math.pow
+import kotlin.math.sqrt
+
+/**
+ * @author tabidachinokaze
+ * @date 2025/11/21
+ */
+/**
+ * 等高线连接器 - 将小线段连接成连续等高线
+ */
+class ContourConnector {
+
+    /**
+     * 连接线段成连续等高线
+     */
+    fun connectSegments(segments: List<List<Point>>): List<List<Point>> {
+        if (segments.isEmpty()) return emptyList()
+
+        val continuousContours = mutableListOf<List<Point>>()
+        val usedSegments = BooleanArray(segments.size) { false }
+
+        for (i in segments.indices) {
+            if (!usedSegments[i] && segments[i].size == 2) {
+                val contour = connectFromSegment(segments, usedSegments, i)
+                if (contour.size >= 2) {
+                    continuousContours.add(contour)
+                }
+            }
+        }
+
+        return continuousContours
+    }
+
+    /**
+     * 从单个线段开始连接
+     */
+    private fun connectFromSegment(
+        segments: List<List<Point>>,
+        usedSegments: BooleanArray,
+        startIndex: Int
+    ): List<Point> {
+        val contour = mutableListOf<Point>()
+        var currentSegment = segments[startIndex]
+
+        // 添加第一个线段
+        contour.addAll(currentSegment)
+        usedSegments[startIndex] = true
+
+        var changed: Boolean
+        do {
+            changed = false
+
+            // 向前连接（在末尾添加）
+            for (i in segments.indices) {
+                if (!usedSegments[i] && segments[i].size == 2) {
+                    val segment = segments[i]
+                    val connection = findConnection(contour.last(), segment)
+                    if (connection != null) {
+                        // 根据连接方向添加点
+                        if (connection == ConnectionType.FORWARD) {
+                            contour.add(segment[1])
+                        } else {
+                            contour.add(segment[0])
+                        }
+                        usedSegments[i] = true
+                        changed = true
+                        break
+                    }
+                }
+            }
+
+            // 向后连接（在开头添加）
+            for (i in segments.indices) {
+                if (!usedSegments[i] && segments[i].size == 2) {
+                    val segment = segments[i]
+                    val connection = findConnection(contour.first(), segment)
+                    if (connection != null) {
+                        // 根据连接方向添加点
+                        if (connection == ConnectionType.FORWARD) {
+                            contour.add(0, segment[0])
+                        } else {
+                            contour.add(0, segment[1])
+                        }
+                        usedSegments[i] = true
+                        changed = true
+                        break
+                    }
+                }
+            }
+
+        } while (changed)
+
+        return contour
+    }
+
+    /**
+     * 查找连接方式
+     */
+    private fun findConnection(point: Point, segment: List<Point>): ConnectionType? {
+        val tolerance = 1e-6 // 连接容差
+
+        if (distance(point, segment[0]) < tolerance) {
+            return ConnectionType.FORWARD
+        }
+        if (distance(point, segment[1]) < tolerance) {
+            return ConnectionType.REVERSE
+        }
+        return null
+    }
+
+    private fun distance(p1: Point, p2: Point): Double {
+        val dx = p1.longitude() - p2.longitude()
+        val dy = p1.latitude() - p2.latitude()
+        return sqrt(dx * dx + dy * dy)
+    }
+
+    private enum class ConnectionType {
+        FORWARD, REVERSE
+    }
+}
+
+/**
+ * 完整的等高线生成器（包含线段连接）
+ */
+class CompleteContourGenerator {
+    private val connector = ContourConnector()
+
+    fun generateContinuousContours(
+        grid: GridModel,
+        contourInterval: Double,
+        minElevation: Double? = null,
+        maxElevation: Double? = null,
+        simplifyTolerance: Double = 0.00001 // 简化容差
+    ): List<ContourLine> {
+        val elevations = grid.cells.filterNotNull()
+        if (elevations.isEmpty()) return emptyList()
+
+        val minElev = minElevation ?: elevations.minOrNull() ?: 0.0
+        val maxElev = maxElevation ?: elevations.maxOrNull() ?: 0.0
+
+        val contours = mutableListOf<ContourLine>()
+
+        var currentElevation = (floor(minElev / contourInterval) + 1) * contourInterval
+        while (currentElevation < maxElev) {
+            // 1. 生成原始线段
+            val rawSegments = marchAllSquares(grid, currentElevation)
+
+            // 2. 连接线段成连续等高线
+            val continuousContours = connector.connectSegments(rawSegments)
+
+            // 3. 简化等高线（可选）
+            val simplifiedContours = continuousContours.map { contour ->
+                simplifyContour(contour, simplifyTolerance)
+            }.filter { it.size >= 2 }
+
+            if (simplifiedContours.isNotEmpty()) {
+                contours.add(ContourLine(currentElevation, simplifiedContours))
+                println("高程 ${currentElevation}m: ${simplifiedContours.size} 条连续等高线")
+            }
+
+            currentElevation += contourInterval
+        }
+
+        return contours
+    }
+
+    /**
+     * 简化等高线（道格拉斯-普克算法）
+     */
+    private fun simplifyContour(contour: List<Point>, tolerance: Double): List<Point> {
+        if (contour.size <= 2) return contour
+
+        return douglasPeucker(contour, 0, contour.size - 1, tolerance)
+    }
+
+    /**
+     * 道格拉斯-普克算法实现
+     */
+    private fun douglasPeucker(
+        points: List<Point>,
+        start: Int,
+        end: Int,
+        tolerance: Double
+    ): List<Point> {
+        if (end <= start + 1) {
+            return listOf(points[start])
+        }
+
+        var maxDistance = 0.0
+        var maxIndex = start
+
+        val startPoint = points[start]
+        val endPoint = points[end]
+
+        // 找到离弦最远的点
+        for (i in start + 1 until end) {
+            val distance = perpendicularDistance(points[i], startPoint, endPoint)
+            if (distance > maxDistance) {
+                maxDistance = distance
+                maxIndex = i
+            }
+        }
+
+        return if (maxDistance > tolerance) {
+            // 递归简化
+            val firstPart = douglasPeucker(points, start, maxIndex, tolerance)
+            val secondPart = douglasPeucker(points, maxIndex, end, tolerance)
+            firstPart.dropLast(1) + secondPart // 避免重复点
+        } else {
+            listOf(startPoint, endPoint)
+        }
+    }
+
+    /**
+     * 计算点到线的垂直距离
+     */
+    private fun perpendicularDistance(point: Point, lineStart: Point, lineEnd: Point): Double {
+        val area = abs(
+            (lineEnd.longitude() - lineStart.longitude()) * (lineStart.latitude() - point.latitude()) -
+                    (lineStart.longitude() - point.longitude()) * (lineEnd.latitude() - lineStart.latitude())
+        )
+        val lineLength = sqrt(
+            (lineEnd.longitude() - lineStart.longitude()).pow(2) +
+                    (lineEnd.latitude() - lineStart.latitude()).pow(2)
+        )
+        return area / lineLength
+    }
+
+    // 原有的移动立方体算法...
+    private fun marchAllSquares(grid: GridModel, elevation: Double): List<List<Point>> {
+        val allLines = mutableListOf<List<Point>>()
+
+        for (r in 0 until grid.rows - 1) {
+            for (c in 0 until grid.cols - 1) {
+                val lines = marchSquare(grid, r, c, elevation)
+                allLines.addAll(lines)
+            }
+        }
+
+        return allLines
+    }
+
+    private fun marchSquare(grid: GridModel, row: Int, col: Int, elevation: Double): List<List<Point>> {
+        val corners = listOf(
+            grid.getValue(row, col),
+            grid.getValue(row, col + 1),
+            grid.getValue(row + 1, col + 1),
+            grid.getValue(row + 1, col)
+        )
+
+        if (corners.any { it == null }) return emptyList()
+
+        val lines = mutableListOf<List<Point>>()
+        val intersections = mutableListOf<Point>()
+
+        // 检查四条边
+        val edges = listOf(
+            Pair(0, 1), // 上
+            Pair(1, 2), // 右
+            Pair(2, 3), // 下
+            Pair(3, 0)  // 左
+        )
+
+        for ((start, end) in edges) {
+            val startElev = corners[start]!!
+            val endElev = corners[end]!!
+
+            if ((startElev - elevation) * (endElev - elevation) < 0) {
+                val t = (elevation - startElev) / (endElev - startElev)
+                val point = calculateEdgePoint(grid, row, col, start, end, t)
+                intersections.add(point)
+            }
+        }
+
+        // 根据交点数量生成线段
+        when (intersections.size) {
+            2 -> lines.add(listOf(intersections[0], intersections[1]))
+            4 -> {
+                // 鞍点情况，需要判断如何连接
+                // 简单处理：按顺序连接
+                lines.add(listOf(intersections[0], intersections[1]))
+                lines.add(listOf(intersections[2], intersections[3]))
+            }
+        }
+
+        return lines
+    }
+
+    private fun calculateEdgePoint(
+        grid: GridModel,
+        row: Int,
+        col: Int,
+        startCorner: Int,
+        endCorner: Int,
+        t: Double
+    ): Point {
+        val cornerPositions = listOf(
+            Pair(col.toDouble(), row.toDouble()),
+            Pair(col + 1.0, row.toDouble()),
+            Pair(col + 1.0, row + 1.0),
+            Pair(col.toDouble(), row + 1.0)
+        )
+
+        val startPos = cornerPositions[startCorner]
+        val endPos = cornerPositions[endCorner]
+
+        val gridX = startPos.first + (endPos.first - startPos.first) * t
+        val gridY = startPos.second + (endPos.second - startPos.second) * t
+
+        return gridToWorld(grid, gridX, gridY)
+    }
+
+    private fun gridToWorld(grid: GridModel, gridX: Double, gridY: Double): Point {
+        val lon = grid.minLon + gridX * (grid.maxLon - grid.minLon) / grid.cols
+        // 修正：Y轴方向取反，因为row=0对应北边（maxLat），row增大向南（minLat）
+        val lat = grid.maxLat - gridY * (grid.maxLat - grid.minLat) / grid.rows
+        return Point.fromLngLat(lon, lat)
+    }
+}
+
+/**
+ * 显示连续等高线
+ */
+fun MapView.displayContinuousContours(
+    grid: GridModel,
+    contourCount: Int = 6,
+    sourceId: String = "continuous-contours",
+    layerId: String = "contour-layer",
+    simplifyTolerance: Double = 0.00001
+) {
+    mapboxMap.getStyle { style ->
+        /*val generator = CompleteContourGenerator()
+        val contours = generator.generateContinuousContours(
+            grid,
+            contourInterval,
+            simplifyTolerance = simplifyTolerance
+        )*/
+
+        val generator = MarchingSquaresContourGenerator()
+        val contours = generator.generateContours(
+            grid,
+            contourCount,
+        )
+
+        val features = mutableListOf<Feature>()
+
+        contours.forEach { contour ->
+            contour.points.forEach { linePoints ->
+                Log.d(
+                    TAG,
+                    buildString {
+                        appendLine("displayContinuousContours: ")
+                        val geoHelper = GeoHelper.getSharedInstance()
+                        val wsg84List = linePoints.map {
+                            Vector3D(x = it.longitude(), y = it.latitude(), z = contour.elevation)
+                        }
+                        appendLine(wsg84List.niceStr())
+                        val points = linePoints.map {
+                            geoHelper.wgs84ToENU(lon = it.longitude(), lat = it.latitude(), hgt = contour.elevation)
+                        }.map {
+                            Vector3D(it.x, it.y, it.z)
+                        }
+                        appendLine(points.niceStr())
+                        appendLine("minX = ${points.minOf { it.x }}")
+                        appendLine("maxX = ${points.maxOf { it.x }}")
+                        appendLine("minY = ${points.minOf { it.y }}")
+                        appendLine("maxY = ${points.maxOf { it.y }}")
+                    }
+                )
+                if (linePoints.size >= 2) {
+                    val lineString = LineString.fromLngLats(linePoints)
+                    val feature = Feature.fromGeometry(lineString)
+
+                    val color = getContourColor(contour.elevation)
+                    feature.addStringProperty("color", color)
+                    feature.addNumberProperty("elevation", contour.elevation.toInt())
+                    feature.addStringProperty("type", "contour-line")
+
+                    features.add(feature)
+                }
+            }
+        }
+
+        setupContourLayers(style, sourceId, layerId, features)
+
+        val totalSegments = contours.sumOf { it.points.size }
+        println("生成 ${contours.size} 条等高线，共 $totalSegments 条连续线段")
+    }
+}

app/src/main/java/com/icegps/geotools/ContourGenerator.kt

@@ -0,0 +1,327 @@
+package com.icegps.geotools
+
+import com.icegps.math.geometry.Vector2D
+import com.mapbox.geojson.Feature
+import com.mapbox.geojson.FeatureCollection
+import com.mapbox.geojson.LineString
+import com.mapbox.geojson.Point
+import com.mapbox.maps.MapView
+import com.mapbox.maps.Style
+import com.mapbox.maps.extension.style.expressions.generated.Expression
+import com.mapbox.maps.extension.style.layers.addLayer
+import com.mapbox.maps.extension.style.layers.generated.LineLayer
+import com.mapbox.maps.extension.style.layers.generated.SymbolLayer
+import com.mapbox.maps.extension.style.layers.properties.generated.SymbolPlacement
+import com.mapbox.maps.extension.style.sources.addSource
+import com.mapbox.maps.extension.style.sources.generated.geoJsonSource
+import kotlin.math.ceil
+
+/**
+ * @author tabidachinokaze
+ * @date 2025/11/21
+ */
+/**
+ * 等高线生成器
+ */
+class ContourGenerator {
+
+    /**
+     * 生成等高线
+     * @param grid 栅格数据
+     * @param contourInterval 等高线间距
+     * @param minElevation 最小高程（可选，自动计算）
+     * @param maxElevation 最大高程（可选，自动计算）
+     */
+    fun generateContours(
+        grid: GridModel,
+        contourInterval: Double,
+        minElevation: Double? = null,
+        maxElevation: Double? = null
+    ): List<ContourLine> {
+        val elevations = grid.cells.filterNotNull()
+        val minElev = minElevation ?: elevations.minOrNull() ?: 0.0
+        val maxElev = maxElevation ?: elevations.maxOrNull() ?: 100.0
+
+        val contours = mutableListOf<ContourLine>()
+
+        // 生成等高线层级
+        var currentElevation = ceil(minElev / contourInterval) * contourInterval
+        while (currentElevation <= maxElev) {
+            val contour = generateContourForElevation(grid, currentElevation)
+            if (contour.isNotEmpty()) {
+                contours.add(ContourLine(currentElevation, contour))
+            }
+            currentElevation += contourInterval
+        }
+
+        return contours
+    }
+
+    /**
+     * 为特定高程生成等高线
+     */
+    private fun generateContourForElevation(
+        grid: GridModel,
+        elevation: Double
+    ): List<List<Point>> {
+        val contours = mutableListOf<List<Point>>()
+        val visited = Array(grid.rows) { BooleanArray(grid.cols) }
+
+        for (r in 0 until grid.rows - 1) {
+            for (c in 0 until grid.cols - 1) {
+                if (!visited[r][c]) {
+                    val contour = traceContour(grid, r, c, elevation, visited)
+                    if (contour.isNotEmpty()) {
+                        contours.add(contour)
+                    }
+                }
+            }
+        }
+
+        return contours
+    }
+
+    /**
+     * 追踪单条等高线
+     */
+    private fun traceContour(
+        grid: GridModel,
+        startRow: Int,
+        startCol: Int,
+        elevation: Double,
+        visited: Array<BooleanArray>
+    ): List<Point> {
+        val contour = mutableListOf<Point>()
+        var row = startRow
+        var col = startCol
+        var direction = 0 // 初始方向
+
+        do {
+            if (row < 0 || row >= grid.rows - 1 || col < 0 || col >= grid.cols - 1) {
+                break
+            }
+
+            if (visited[row][col]) {
+                break
+            }
+
+            visited[row][col] = true
+
+            // 获取当前网格的四个角点高程
+            val corners = listOf(
+                grid.getValue(row, col) ?: continue,
+                grid.getValue(row, col + 1) ?: continue,
+                grid.getValue(row + 1, col + 1) ?: continue,
+                grid.getValue(row + 1, col) ?: continue
+            )
+
+            // 计算网格边上的交点
+            val intersections = calculateVector2Ds(grid, row, col, elevation, corners)
+
+            if (intersections.isNotEmpty()) {
+                // 选择第一个交点（简化处理）
+                val intersection = intersections.first()
+                val point = gridToWorld(grid, intersection.x, intersection.y)
+                contour.add(point)
+
+                // 移动到下一个网格（简化版，实际应该根据交点位置决定方向）
+                when (direction) {
+                    0 -> col++ // 右
+                    1 -> row++ // 下
+                    2 -> col-- // 左  
+                    3 -> row-- // 上
+                }
+                direction = (direction + 1) % 4
+            } else {
+                break
+            }
+
+        } while (row != startRow || col != startCol)
+
+        return contour
+    }
+
+    /**
+     * 计算网格边上的交点
+     */
+    private fun calculateVector2Ds(
+        grid: GridModel,
+        row: Int,
+        col: Int,
+        elevation: Double,
+        corners: List<Double>
+    ): List<Vector2D> {
+        val intersections = mutableListOf<Vector2D>()
+        val (a, b, c, d) = corners
+
+        // 检查四条边
+        // 上边 (a-b)
+        if ((a - elevation) * (b - elevation) < 0) {
+            val t = (elevation - a) / (b - a)
+            intersections.add(Vector2D(col + t, row.toDouble()))
+        }
+
+        // 右边 (b-c)
+        if ((b - elevation) * (c - elevation) < 0) {
+            val t = (elevation - b) / (c - b)
+            intersections.add(Vector2D(col + 1.0, row + t))
+        }
+
+        // 下边 (c-d)
+        if ((c - elevation) * (d - elevation) < 0) {
+            val t = (elevation - c) / (d - c)
+            intersections.add(Vector2D(col + 1 - t, row + 1.0))
+        }
+
+        // 左边 (d-a)
+        if ((d - elevation) * (a - elevation) < 0) {
+            val t = (elevation - d) / (a - d)
+            intersections.add(Vector2D(col.toDouble(), row + 1 - t))
+        }
+
+        return intersections
+    }
+
+    /**
+     * 网格坐标转世界坐标
+     */
+    private fun gridToWorld(grid: GridModel, gridX: Double, gridY: Double): Point {
+        val lon = grid.minLon + gridX * (grid.maxLon - grid.minLon) / grid.cols
+        // 修正：Y轴方向取反，因为row=0对应北边（maxLat），row增大向南（minLat）
+        val lat = grid.maxLat - gridY * (grid.maxLat - grid.minLat) / grid.rows
+        return Point.fromLngLat(lon, lat)
+    }
+}
+
+/**
+ * 等高线数据类
+ */
+data class ContourLine(
+    val elevation: Double,          // 等高线高程
+    val points: List<List<Point>>,  // 等高线路径（可能有多个闭合环）
+    val properties: Map<String, Any> = emptyMap()
+) {
+    fun getLineStrings(): List<LineString> {
+        return points.map { LineString.fromLngLats(it) }
+    }
+}
+
+// ========================================绘制等高线==========================================
+/**
+ * 在Mapbox中显示等高线
+ */
+fun MapView.displayContours(
+    grid: GridModel,
+    contourInterval: Double = 10.0,
+    sourceId: String = "contour-lines",
+    layerId: String = "contour-layer"
+) {
+    mapboxMap.getStyle { style ->
+        // val generator = ContourGenerator()
+        // val contours = generator.generateContours(grid, contourInterval)
+        val simpleContourGenerator = SimpleContourGenerator()
+        val contours = simpleContourGenerator.generateSimpleContours(
+            grid = grid,
+            contourInterval = contourInterval
+        )
+
+        val features = mutableListOf<Feature>()
+
+        // 为每条等高线创建要素
+        contours.forEach { contour ->
+            contour.points.forEach { linePoints ->
+                if (linePoints.size >= 2) {
+                    val lineString = LineString.fromLngLats(linePoints)
+                    val feature = Feature.fromGeometry(lineString)
+
+                    // 根据高程设置颜色和属性
+                    val color = getContourColor(contour.elevation)
+                    feature.addStringProperty("color", color)
+                    feature.addNumberProperty("elevation", contour.elevation)
+                    feature.addStringProperty("type", "contour-line")
+
+                    features.add(feature)
+                }
+            }
+        }
+
+        setupContourLayers(style, sourceId, layerId, features)
+
+        println("生成 ${contours.size} 条等高线，共 ${features.size} 条线段")
+    }
+}
+
+/**
+ * 根据高程获取等高线颜色
+ */
+fun getContourColor(elevation: Double): String {
+    return when {
+        elevation < 100 -> "#FF0000" // 红色 - 低海拔
+        elevation < 200 -> "#FF9900" // 橙色
+        elevation < 300 -> "#FFFF00" // 黄色
+        elevation < 400 -> "#00FF00" // 绿色
+        elevation < 500 -> "#00FFFF" // 青色
+        else -> "#0000FF"           // 蓝色 - 高海拔
+    }
+}
+
+/**
+ * 设置等高线图层
+ */
+fun MapView.setupContourLayers(
+    style: Style,
+    sourceId: String,
+    layerId: String,
+    features: List<Feature>
+) {
+    val source = geoJsonSource(sourceId) {
+        featureCollection(FeatureCollection.fromFeatures(features))
+    }
+
+    // 清理旧图层
+    try {
+        style.removeStyleLayer(layerId)
+    } catch (_: Exception) {
+    }
+    try {
+        style.removeStyleLayer("$layerId-labels")
+    } catch (_: Exception) {
+    }
+
+    if (style.styleSourceExists(sourceId)) {
+        style.removeStyleSource(sourceId)
+    }
+
+    // 添加数据源
+    style.addSource(source)
+
+    // 等高线图层
+    val lineLayer = LineLayer(layerId, sourceId).apply {
+        lineColor(Expression.toColor(Expression.get("color")))
+        lineWidth(
+            Expression.interpolate(
+                Expression.exponential(1.0),
+                Expression.zoom(),
+                Expression.literal(10.0), Expression.literal(0.5),  // 缩放级别10，线宽0.5
+                Expression.literal(15.0), Expression.literal(1.0),  // 缩放级别15，线宽1.0
+                Expression.literal(20.0), Expression.literal(2.0)   // 缩放级别20，线宽2.0
+            )
+        )
+        lineOpacity(0.8)
+        filter(Expression.eq(Expression.get("type"), Expression.literal("contour-line")))
+    }
+    style.addLayer(lineLayer)
+
+    // 等高线标注图层（可选）
+    val labelLayer = SymbolLayer("$layerId-labels", sourceId).apply {
+        textField(Expression.toString(Expression.get("elevation")))
+        textSize(10.0)
+        textColor("#000000")
+        textHaloColor("#FFFFFF")
+        textHaloWidth(1.0)
+        textOpacity(0.9)
+        symbolPlacement(SymbolPlacement.LINE)
+        filter(Expression.eq(Expression.get("type"), Expression.literal("contour-line")))
+    }
+    style.addLayer(labelLayer)
+}

app/src/main/java/com/icegps/geotools/DraggableTrendArrow.kt

@@ -0,0 +1,314 @@
+package com.icegps.geotools
+
+import com.mapbox.android.gestures.MoveGestureDetector
+import com.mapbox.geojson.Feature
+import com.mapbox.geojson.LineString
+import com.mapbox.geojson.Point
+import com.mapbox.geojson.Polygon
+import com.mapbox.maps.MapView
+import com.mapbox.maps.plugin.gestures.OnMoveListener
+import com.mapbox.maps.plugin.gestures.addOnMoveListener
+import com.mapbox.maps.plugin.gestures.gestures
+import com.mapbox.maps.plugin.gestures.removeOnMoveListener
+import kotlin.math.*
+
+/**
+ * @author tabidachinokaze
+ * @date 2025/11/20
+ */
+/**
+ * 可拖动的整体趋势箭头
+ */
+class DraggableTrendArrow(
+    private val mapView: MapView,
+    private val grid: GridModel,
+    private val initialSlopeDirection: Double
+) {
+    private var currentDirection = initialSlopeDirection
+    private val sourceId = "draggable-trend-arrow"
+    private val lineLayerId = "draggable-trend-line"
+    private val headLayerId = "draggable-trend-head"
+
+    private val centerLon = (grid.minLon + grid.maxLon) / 2
+    private val centerLat = (grid.minLat + grid.maxLat) / 2
+    private val baseArrowLength = (grid.maxLon - grid.minLon) * 0.3
+
+    // 回调函数
+    var onDirectionChanged: ((Double) -> Unit)? = null
+
+    /**
+     * 显示可拖动的箭头
+     */
+    fun show() {
+        updateArrow()
+        setupGestureListeners()
+    }
+
+    /**
+     * 隐藏箭头
+     */
+    fun hide() {
+        mapView.mapboxMap.getStyle { style ->
+            try {
+                style.removeStyleLayer(lineLayerId)
+            } catch (_: Exception) {
+            }
+            try {
+                style.removeStyleLayer(headLayerId)
+            } catch (_: Exception) {
+            }
+            if (style.styleSourceExists(sourceId)) {
+                style.removeStyleSource(sourceId)
+            }
+        }
+    }
+
+    /**
+     * 更新箭头方向
+     */
+    fun updateDirection(newDirection: Double) {
+        currentDirection = newDirection
+        updateArrow()
+        onDirectionChanged?.invoke(newDirection)
+    }
+
+    /**
+     * 获取当前方向
+     */
+    fun getCurrentDirection(): Double = currentDirection
+
+    /**
+     * 更新箭头显示
+     */
+    private fun updateArrow() {
+        mapView.mapboxMap.getStyle { style ->
+            val arrowDirectionRad = Math.toRadians(currentDirection)
+            val endLon = centerLon + cos(arrowDirectionRad) * baseArrowLength
+            val endLat = centerLat + sin(arrowDirectionRad) * baseArrowLength
+
+            // 创建箭杆
+            val arrowLine = LineString.fromLngLats(
+                listOf(
+                    Point.fromLngLat(centerLon, centerLat),
+                    Point.fromLngLat(endLon, endLat)
+                )
+            )
+
+            // 创建箭头头部
+            val headSize = baseArrowLength * 0.15
+            val leftRad = arrowDirectionRad + Math.PI * 0.8
+            val rightRad = arrowDirectionRad - Math.PI * 0.8
+
+            val leftLon = endLon + cos(leftRad) * headSize
+            val leftLat = endLat + sin(leftRad) * headSize
+            val rightLon = endLon + cos(rightRad) * headSize
+            val rightLat = endLat + sin(rightRad) * headSize
+
+            val headRing = listOf(
+                Point.fromLngLat(endLon, endLat),
+                Point.fromLngLat(leftLon, leftLat),
+                Point.fromLngLat(rightLon, rightLat),
+                Point.fromLngLat(endLon, endLat)
+            )
+            val headPolygon = Polygon.fromLngLats(listOf(headRing))
+
+            val features = listOf(
+                Feature.fromGeometry(arrowLine).apply {
+                    addStringProperty("color", "#FF6B35")
+                    addStringProperty("type", "draggable-arrow-line")
+                },
+                Feature.fromGeometry(headPolygon).apply {
+                    addStringProperty("color", "#FF6B35")
+                    addStringProperty("type", "draggable-arrow-head")
+                }
+            )
+
+            // 更新或创建图层
+            // updateOrCreateArrowLayer(style, features)
+
+            println("更新箭头方向: ${"%.1f".format(currentDirection)}°")
+        }
+    }
+
+    /**
+     * 设置手势监听器
+     */
+    private fun setupGestureListeners() {
+        // 添加点击检测
+        mapView.gestures.addOnMapClickListener { point ->
+            if (isPointOnArrow(point)) {
+                startDragging(point)
+                true
+            } else {
+                false
+            }
+        }
+
+        // 添加长按检测（更易触发）
+        mapView.gestures.addOnMapLongClickListener { point ->
+            if (isPointOnArrow(point)) {
+                startDragging(point)
+                true
+            } else {
+                false
+            }
+        }
+    }
+
+    /**
+     * 检查点击点是否在箭头上
+     */
+    private fun isPointOnArrow(point: Point): Boolean {
+        val arrowDirectionRad = Math.toRadians(currentDirection)
+        val endLon = centerLon + cos(arrowDirectionRad) * baseArrowLength
+        val endLat = centerLat + sin(arrowDirectionRad) * baseArrowLength
+
+        // 计算点击点与箭头线的距离
+        val distance = calculateDistanceToLine(
+            point.longitude(), point.latitude(),
+            centerLon, centerLat, endLon, endLat
+        )
+
+        // 如果距离小于阈值，认为点击在箭头上
+        val threshold = 0.0002 // 约20米
+        return distance < threshold
+    }
+
+    /**
+     * 开始拖动
+     */
+    private fun startDragging(startPoint: Point) {
+        var isDragging = true
+        val moveListener = object : OnMoveListener {
+            override fun onMove(detector: MoveGestureDetector): Boolean {
+                return true
+            }
+
+            override fun onMoveBegin(detector: MoveGestureDetector) {
+                if (isDragging) {
+                    val screenPoint = detector.focalPoint
+                    // val point = detector.point
+                    // updateArrowDirectionFromPoint(point)
+                }
+            }
+
+            override fun onMoveEnd(detector: MoveGestureDetector) {
+                isDragging = false
+                mapView.mapboxMap.removeOnMoveListener(this)
+            }
+        }
+        // 触摸移动监听
+        mapView.mapboxMap.addOnMoveListener(moveListener)
+    }
+
+    /**
+     * 根据拖动点更新箭头方向
+     */
+    private fun updateArrowDirectionFromPoint(point: Point) {
+        val dx = point.longitude() - centerLon
+        val dy = point.latitude() - centerLat
+
+        // 计算新方向（弧度）
+        var newDirectionRad = atan2(dy, dx)
+
+        // 转换为角度并确保在 0-360 范围内
+        var newDirection = Math.toDegrees(newDirectionRad)
+        if (newDirection < 0) newDirection += 360.0
+
+        updateDirection(newDirection)
+    }
+
+    /**
+     * 计算点到线的距离
+     */
+    private fun calculateDistanceToLine(
+        px: Double, py: Double,
+        x1: Double, y1: Double,
+        x2: Double, y2: Double
+    ): Double {
+        val A = px - x1
+        val B = py - y1
+        val C = x2 - x1
+        val D = y2 - y1
+
+        val dot = A * C + B * D
+        val lenSq = C * C + D * D
+        var param = -1.0
+
+        if (lenSq != 0.0) {
+            param = dot / lenSq
+        }
+
+        var xx: Double
+        var yy: Double
+
+        if (param < 0) {
+            xx = x1
+            yy = y1
+        } else if (param > 1) {
+            xx = x2
+            yy = y2
+        } else {
+            xx = x1 + param * C
+            yy = y1 + param * D
+        }
+
+        val dx = px - xx
+        val dy = py - yy
+        return sqrt(dx * dx + dy * dy)
+    }
+
+    /**
+     * 更新或创建箭头图层
+     */
+    /*private fun updateOrCreateArrowLayer(style: Style, features: List<Feature>) {
+        val source = geoJsonSource(sourceId) {
+            featureCollection(FeatureCollection.fromFeatures(features))
+        }
+
+        if (style.styleSourceExists(sourceId)) {
+            // 更新现有数据源
+            (style.getStyleSource(sourceId) as GeoJsonSource).featureCollection(
+                featureCollection(FeatureCollection.fromFeatures(features))
+            )
+        } else {
+            // 创建新数据源和图层
+            style.addSource(source)
+
+            // 箭杆图层
+            val lineLayer = LineLayer(lineLayerId, sourceId).apply {
+                lineColor(Expression.toColor(Expression.get("color")))
+                lineWidth(4.0)
+                lineCap(LineCap.ROUND)
+                withFilter(Expression.eq(Expression.get("type"), Expression.literal("draggable-arrow-line")))
+            }
+            style.addLayer(lineLayer)
+
+            // 箭头头部图层
+            val headLayer = FillLayer(headLayerId, sourceId).apply {
+                fillColor(Expression.toColor(Expression.get("color")))
+                fillOpacity(1.0)
+                withFilter(Expression.eq(Expression.get("type"), Expression.literal("draggable-arrow-head")))
+            }
+            style.addLayer(headLayer)
+        }
+    }*/
+}
+
+/**
+ * 简化的可拖动箭头显示函数
+ */
+fun MapView.displayDraggableTrendArrow(
+    grid: GridModel,
+    initialDirection: Double? = null,
+    onDirectionChanged: (Double) -> Unit = {}
+): DraggableTrendArrow {
+    // 如果没有提供初始方向，计算整体趋势
+    val direction = initialDirection ?: grid.calculateOverallSlopeTrend().direction
+
+    val draggableArrow = DraggableTrendArrow(this, grid, direction)
+    draggableArrow.onDirectionChanged = onDirectionChanged
+    draggableArrow.show()
+
+    return draggableArrow
+}

app/src/main/java/com/icegps/geotools/EarthworkManager.kt

@@ -0,0 +1,937 @@
+package com.icegps.geotools
+
+import android.util.Log
+import com.mapbox.geojson.Feature
+import com.mapbox.geojson.FeatureCollection
+import com.mapbox.geojson.LineString
+import com.mapbox.geojson.Point
+import com.mapbox.geojson.Polygon
+import com.mapbox.maps.MapView
+import com.mapbox.maps.Style
+import com.mapbox.maps.extension.style.expressions.generated.Expression
+import com.mapbox.maps.extension.style.layers.addLayer
+import com.mapbox.maps.extension.style.layers.generated.FillLayer
+import com.mapbox.maps.extension.style.layers.generated.LineLayer
+import com.mapbox.maps.extension.style.layers.generated.SymbolLayer
+import com.mapbox.maps.extension.style.layers.properties.generated.LineCap
+import com.mapbox.maps.extension.style.layers.properties.generated.LineJoin
+import com.mapbox.maps.extension.style.sources.addSource
+import com.mapbox.maps.extension.style.sources.generated.geoJsonSource
+import kotlin.math.abs
+import kotlin.math.cos
+import kotlin.math.sin
+import kotlin.math.sqrt
+
+/**
+ * @author tabidachinokaze
+ * @date 2025/11/20
+ */
+/**
+ * 填平计算器 - 将区域填平到平均高度
+ */
+class LevelingCalculator {
+
+    /**
+     * 计算填平到平均高度的土方量
+     * 平均高度 = 所有点高程的平均值
+     */
+    fun calculateLeveling(grid: GridModel): LevelingResult {
+        val elevations = grid.cells.filterNotNull()
+        if (elevations.isEmpty()) {
+            throw IllegalArgumentException("网格中没有有效的高程数据")
+        }
+
+        val averageElevation = elevations.average()
+        val calculator = EarthworkCalculator()
+        val earthworkResult = calculator.calculateForFlatDesign(grid, averageElevation)
+
+        return LevelingResult(
+            targetElevation = averageElevation,
+            earthworkResult = earthworkResult,
+            elevationStats = ElevationStats(
+                min = elevations.minOrNull() ?: 0.0,
+                max = elevations.maxOrNull() ?: 0.0,
+                mean = averageElevation,
+                stdDev = calculateStdDev(elevations),
+                validCount = elevations.size
+            )
+        )
+    }
+
+    /**
+     * 计算填平到指定高度的土方量
+     */
+    fun calculateLevelingToElevation(grid: GridModel, targetElevation: Double): LevelingResult {
+        val elevations = grid.cells.filterNotNull()
+        val calculator = EarthworkCalculator()
+        val earthworkResult = calculator.calculateForFlatDesign(grid, targetElevation)
+
+        return LevelingResult(
+            targetElevation = targetElevation,
+            earthworkResult = earthworkResult,
+            elevationStats = ElevationStats(
+                min = elevations.minOrNull() ?: 0.0,
+                max = elevations.maxOrNull() ?: 0.0,
+                mean = elevations.average(),
+                stdDev = calculateStdDev(elevations),
+                validCount = elevations.size
+            )
+        )
+    }
+}
+
+data class LevelingResult(
+    val targetElevation: Double,      // 目标高程
+    val earthworkResult: EarthworkResult, // 土方量结果
+    val elevationStats: ElevationStats    // 原始高程统计
+) {
+    override fun toString(): String {
+        return "填平到 ${"%.2f".format(targetElevation)} m:\n" +
+                "原始高程: ${"%.1f".format(elevationStats.min)}-${"%.1f".format(elevationStats.max)} m " +
+                "(平均${"%.1f".format(elevationStats.mean)} m)\n" +
+                "挖方: ${"%.0f".format(earthworkResult.cutVolume)} m³, " +
+                "填方: ${"%.0f".format(earthworkResult.fillVolume)} m³, " +
+                "净土方: ${"%.0f".format(earthworkResult.netVolume)} m³"
+    }
+}
+
+/**
+ * 斜坡计算器 - 在区域内创建指定坡向的斜面
+ */
+class SlopeCalculator {
+
+    /**
+     * 计算斜坡方案的土方量
+     * @param grid 原始地形网格
+     * @param slopeDirection 坡向 (度，0=北，90=东)
+     * @param slopePercentage 坡度 (%)
+     * @param baseHeightOffset 基准面高度偏移 (m)，正数上移，负数下移
+     */
+    fun calculateSlope(
+        grid: GridModel,
+        slopeDirection: Double,
+        slopePercentage: Double,
+        baseHeightOffset: Double = 0.0
+    ): SlopeResult {
+        // 计算区域的中心点作为基准点
+        val centerLon = (grid.minLon + grid.maxLon) / 2
+        val centerLat = (grid.minLat + grid.maxLat) / 2
+
+        // 计算基准点的高程（使用平均高程 + 偏移）
+        val elevations = grid.cells.filterNotNull()
+        val baseElevation = elevations.average() + baseHeightOffset
+
+        val basePoint = Triple(centerLon, centerLat, baseElevation)
+        val calculator = EarthworkCalculator()
+        val earthworkResult = calculator.calculateForSlopeDesign(
+            grid, basePoint, slopePercentage, slopeDirection
+        )
+
+        return SlopeResult(
+            slopeDirection = slopeDirection,
+            slopePercentage = slopePercentage,
+            baseHeightOffset = baseHeightOffset,
+            baseElevation = baseElevation,
+            earthworkResult = earthworkResult,
+            designSurface = generateSlopeDesignGrid(grid, basePoint, slopePercentage, slopeDirection)
+        )
+    }
+
+    /**
+     * 生成斜坡设计面网格（用于可视化）
+     */
+    private fun generateSlopeDesignGrid(
+        grid: GridModel,
+        basePoint: Triple<Double, Double, Double>,
+        slopePercentage: Double,
+        slopeDirection: Double
+    ): GridModel {
+        val designCells = Array<Double?>(grid.rows * grid.cols) { null }
+        val (baseLon, baseLat, baseElev) = basePoint
+        val slopeRatio = slopePercentage / 100.0
+
+        for (r in 0 until grid.rows) {
+            for (c in 0 until grid.cols) {
+                if (grid.getValue(r, c) != null) {
+                    val cellLon = grid.minLon + (c + 0.5) * (grid.maxLon - grid.minLon) / grid.cols
+                    val cellLat = grid.minLat + (r + 0.5) * (grid.maxLat - grid.minLat) / grid.rows
+
+                    val designElev = calculateSlopeElevation(
+                        cellLon, cellLat, baseLon, baseLat, baseElev, slopeRatio, slopeDirection
+                    )
+                    designCells[r * grid.cols + c] = designElev
+                }
+            }
+        }
+
+        return GridModel(
+            minLon = grid.minLon,
+            minLat = grid.minLat,
+            maxLon = grid.maxLon,
+            maxLat = grid.maxLat,
+            rows = grid.rows,
+            cols = grid.cols,
+            cellSizeMeters = grid.cellSizeMeters,
+            cells = designCells
+        )
+    }
+
+    /**
+     * 计算斜坡上某点的设计高程
+     */
+    private fun calculateSlopeElevationLegecy(
+        pointLon: Double, pointLat: Double,
+        baseLon: Double, baseLat: Double, baseElev: Double,
+        slopeRatio: Double, aspect: Double
+    ): Double {
+        val dx = (pointLon - baseLon) * 111320.0 * cos(Math.toRadians(baseLat))
+        val dy = (pointLat - baseLat) * 111320.0
+        val aspectRad = Math.toRadians(aspect)
+        val distanceInSlopeDirection = dx * cos(aspectRad) + dy * sin(aspectRad)
+        val heightDiff = distanceInSlopeDirection * slopeRatio
+        return baseElev + heightDiff
+    }
+
+    /**
+     * 修正的斜坡高程计算
+     */
+    private fun calculateSlopeElevation(
+        pointLon: Double, pointLat: Double,
+        baseLon: Double, baseLat: Double, baseElev: Double,
+        slopeRatio: Double, slopeDirection: Double
+    ): Double {
+        // 计算相对位置向量（在ENU坐标系中，Y指向南）
+        val east = (pointLon - baseLon) * 111320.0 * cos(Math.toRadians(baseLat))
+        val south = (pointLat - baseLat) * 111320.0  // 注意：这是南方向
+
+        val slopeRad = Math.toRadians(slopeDirection)
+        val slopeVectorX = cos(slopeRad)  // 东方向分量
+        val slopeVectorY = -sin(slopeRad) // 南方向分量（注意负号）
+
+        // 计算在坡度方向上的投影距离
+        val projection = east * slopeVectorX + south * slopeVectorY
+        val heightDiff = projection * slopeRatio
+
+        return baseElev + heightDiff
+    }
+}
+
+data class SlopeResult(
+    val slopeDirection: Double,           // 坡向 (度)
+    val slopePercentage: Double,          // 坡度 (%)
+    val baseHeightOffset: Double,         // 基准面高度偏移 (m)
+    val baseElevation: Double,            // 基准点高程 (m)
+    val earthworkResult: EarthworkResult, // 土方量结果
+    val designSurface: GridModel          // 设计面网格（用于可视化）
+) {
+    fun getDirectionName(): String {
+        return when {
+            slopeDirection >= 337.5 || slopeDirection < 22.5 -> "北"
+            slopeDirection >= 22.5 && slopeDirection < 67.5 -> "东北"
+            slopeDirection >= 67.5 && slopeDirection < 112.5 -> "东"
+            slopeDirection >= 112.5 && slopeDirection < 157.5 -> "东南"
+            slopeDirection >= 157.5 && slopeDirection < 202.5 -> "南"
+            slopeDirection >= 202.5 && slopeDirection < 247.5 -> "西南"
+            slopeDirection >= 247.5 && slopeDirection < 292.5 -> "西"
+            else -> "西北"
+        }
+    }
+
+    override fun toString(): String {
+        return "斜坡设计: ${getDirectionName()}方向 ${"%.1f".format(slopePercentage)}% 坡度\n" +
+                "基准高程: ${"%.2f".format(baseElevation)} m (偏移 ${"%.2f".format(baseHeightOffset)} m)\n" +
+                "挖方: ${"%.0f".format(earthworkResult.cutVolume)} m³, " +
+                "填方: ${"%.0f".format(earthworkResult.fillVolume)} m³, " +
+                "净土方: ${"%.0f".format(earthworkResult.netVolume)} m³"
+    }
+}
+
+/**
+ * 土方量计算管理器
+ */
+class EarthworkManager {
+
+    private val levelingCalculator = LevelingCalculator()
+    private val slopeCalculator = SlopeCalculator()
+
+    /**
+     * 需求1: 填平计算
+     */
+    fun calculateLeveling(grid: GridModel, targetElevation: Double? = null): LevelingResult {
+        return if (targetElevation != null) {
+            levelingCalculator.calculateLevelingToElevation(grid, targetElevation)
+        } else {
+            levelingCalculator.calculateLeveling(grid)
+        }
+    }
+
+    /**
+     * 需求2: 斜坡计算
+     */
+    fun calculateSlope(
+        grid: GridModel,
+        slopeDirection: Double,
+        slopePercentage: Double = 2.0,
+        baseHeightOffset: Double = 0.0
+    ): SlopeResult {
+        return slopeCalculator.calculateSlope(grid, slopeDirection, slopePercentage, baseHeightOffset)
+    }
+
+    /**
+     * 比较多种方案
+     */
+    fun compareSolutions(grid: GridModel, slopeDirection: Double): List<SolutionComparison> {
+        val solutions = mutableListOf<SolutionComparison>()
+
+        // 方案1: 填平到平均高度
+        val levelingResult = calculateLeveling(grid)
+        solutions.add(
+            SolutionComparison(
+                type = "填平到平均高度",
+                description = "将整个区域填平到平均高程",
+                result = levelingResult.earthworkResult,
+                additionalInfo = "目标高程: ${"%.2f".format(levelingResult.targetElevation)} m"
+            )
+        )
+
+        // 方案2: 斜坡方案（当前坡向）
+        val slopeResult = calculateSlope(grid, slopeDirection)
+        solutions.add(
+            SolutionComparison(
+                type = "斜坡设计",
+                description = "${slopeResult.getDirectionName()}方向 ${"%.1f".format(slopeResult.slopePercentage)}% 坡度",
+                result = slopeResult.earthworkResult,
+                additionalInfo = "基准高程: ${"%.2f".format(slopeResult.baseElevation)} m"
+            )
+        )
+
+        // 方案3: 填平到最优高度（净土方量最小）
+        val optimalResult = findOptimalLeveling(grid)
+        solutions.add(
+            SolutionComparison(
+                type = "填平到最优高度",
+                description = "使净土方量最小的填平高度",
+                result = optimalResult.earthworkResult,
+                additionalInfo = "最优高程: ${"%.2f".format(optimalResult.targetElevation)} m"
+            )
+        )
+
+        return solutions.sortedBy { abs(it.result.netVolume) }
+    }
+
+    /**
+     * 寻找使净土方量最小的填平高度
+     */
+    private fun findOptimalLeveling(grid: GridModel): LevelingResult {
+        val elevations = grid.cells.filterNotNull()
+        val minElev = elevations.minOrNull() ?: 0.0
+        val maxElev = elevations.maxOrNull() ?: 100.0
+
+        var bestElevation = minElev
+        var bestNetVolume = Double.MAX_VALUE
+
+        // 在最小和最大高程之间搜索
+        for (elev in minElev.toInt()..maxElev.toInt()) {
+            val result = levelingCalculator.calculateLevelingToElevation(grid, elev.toDouble())
+            val netVolume = abs(result.earthworkResult.netVolume)
+
+            if (netVolume < bestNetVolume) {
+                bestNetVolume = netVolume
+                bestElevation = elev.toDouble()
+            }
+        }
+
+        return levelingCalculator.calculateLevelingToElevation(grid, bestElevation)
+    }
+}
+
+data class SolutionComparison(
+    val type: String,
+    val description: String,
+    val result: EarthworkResult,
+    val additionalInfo: String
+)
+
+/**
+ * 土方量计算结果
+ */
+data class EarthworkResult(
+    val cutVolume: Double,    // 挖方量 (m³)
+    val fillVolume: Double,   // 填方量 (m³)
+    val netVolume: Double,    // 净土方量 (m³)
+    val cutArea: Double,      // 挖方面积 (m²)
+    val fillArea: Double,     // 填方面积 (m²)
+    val totalArea: Double     // 总面积 (m²)
+) {
+    override fun toString(): String {
+        return buildString {
+            appendLine("EarthworkResult")
+            appendLine("挖方: ${"%.1f".format(cutVolume)} m³")
+            appendLine("填方: ${"%.1f".format(fillVolume)} m³")
+            appendLine("净土方: ${"%.1f".format(netVolume)} m³")
+            appendLine("挖方面积: ${"%.1f".format(cutArea)} m²")
+            appendLine("填方面积: ${"%.1f".format(fillArea)} m²")
+        }
+    }
+}
+
+/**
+ * 基于GridModel的土方量计算器
+ */
+class EarthworkCalculator {
+
+    /**
+     * 计算平面设计方案的土方量
+     * @param designElevation 设计高程
+     */
+    fun calculateForFlatDesign(
+        grid: GridModel,
+        designElevation: Double
+    ): EarthworkResult {
+        var cutVolume = 0.0
+        var fillVolume = 0.0
+        var cutArea = 0.0
+        var fillArea = 0.0
+        val cellArea = grid.cellSizeMeters * grid.cellSizeMeters
+
+        for (r in 0 until grid.rows) {
+            for (c in 0 until grid.cols) {
+                val originalElev = grid.getValue(r, c) ?: continue
+
+                val heightDiff = designElevation - originalElev
+                val volume = heightDiff * cellArea
+
+                if (volume > 0) {
+                    fillVolume += volume
+                    fillArea += cellArea
+                } else if (volume < 0) {
+                    cutVolume += abs(volume)
+                    cutArea += cellArea
+                }
+                // volume == 0 时不计算面积
+            }
+        }
+
+        val totalArea = (cutArea + fillArea) / 10000.0 // 转换为公顷
+
+        return EarthworkResult(
+            cutVolume = cutVolume,
+            fillVolume = fillVolume,
+            netVolume = fillVolume - cutVolume,
+            cutArea = cutArea,
+            fillArea = fillArea,
+            totalArea = totalArea
+        )
+    }
+
+    /**
+     * 计算斜面设计方案的土方量
+     * @param basePoint 基准点 (lon, lat, elevation)
+     * @param slope 坡度 (%)
+     * @param aspect 坡向 (度，0=北，90=东)
+     */
+    fun calculateForSlopeDesign(
+        grid: GridModel,
+        basePoint: Triple<Double, Double, Double>,
+        slope: Double,
+        aspect: Double
+    ): EarthworkResult {
+        var cutVolume = 0.0
+        var fillVolume = 0.0
+        var cutArea = 0.0
+        var fillArea = 0.0
+        val cellArea = grid.cellSizeMeters * grid.cellSizeMeters
+
+        val (baseLon, baseLat, baseElev) = basePoint
+        val slopeRatio = slope / 100.0 // 转换为小数
+
+        for (r in 0 until grid.rows) {
+            for (c in 0 until grid.cols) {
+                val originalElev = grid.getValue(r, c) ?: continue
+
+                // 计算当前网格中心坐标
+                val cellLon = grid.minLon + (c + 0.5) * (grid.maxLon - grid.minLon) / grid.cols
+                val cellLat = grid.minLat + (r + 0.5) * (grid.maxLat - grid.minLat) / grid.rows
+
+                // 计算设计高程
+                val designElev = calculateSlopeElevation(
+                    cellLon, cellLat, baseLon, baseLat, baseElev, slopeRatio, aspect
+                )
+
+                val heightDiff = designElev - originalElev
+                val volume = heightDiff * cellArea
+
+                if (volume > 0) {
+                    fillVolume += volume
+                    fillArea += cellArea
+                } else if (volume < 0) {
+                    cutVolume += abs(volume)
+                    cutArea += cellArea
+                }
+            }
+        }
+
+        val totalArea = (cutArea + fillArea) / 10000.0
+
+        return EarthworkResult(
+            cutVolume = cutVolume,
+            fillVolume = fillVolume,
+            netVolume = fillVolume - cutVolume,
+            cutArea = cutArea,
+            fillArea = fillArea,
+            totalArea = totalArea
+        )
+    }
+
+    /**
+     * 计算斜面某点的设计高程
+     */
+    private fun calculateSlopeElevation(
+        pointLon: Double, pointLat: Double,
+        baseLon: Double, baseLat: Double, baseElev: Double,
+        slopeRatio: Double, aspect: Double
+    ): Double {
+        // 计算点到基准点的水平距离（近似计算）
+        val dx = (pointLon - baseLon) * 111320.0 * cos(Math.toRadians(baseLat)) // 东西距离
+        val dy = (pointLat - baseLat) * 111320.0 // 南北距离
+
+        // 计算在坡向方向上的投影距离
+        val aspectRad = Math.toRadians(aspect)
+        val distanceInSlopeDirection = dx * cos(aspectRad) + dy * sin(aspectRad)
+
+        // 计算高差
+        val heightDiff = distanceInSlopeDirection * slopeRatio
+
+        return baseElev + heightDiff
+    }
+}
+
+data class ElevationStats(
+    val min: Double,
+    val max: Double,
+    val mean: Double,
+    val stdDev: Double,
+    val validCount: Int
+)
+
+private fun calculateStdDev(values: List<Double>): Double {
+    if (values.isEmpty()) return 0.0
+    val mean = values.average()
+    val variance = values.map { (it - mean) * (it - mean) }.average()
+    return sqrt(variance)
+}
+
+/*====================================================================================*/
+
+/**
+ * 在Mapbox上绘制填平结果
+ */
+fun MapView.displayLevelingResult(
+    originalGrid: GridModel,
+    levelingResult: LevelingResult,
+    sourceId: String = "leveling-result",
+    layerId: String = "leveling-layer"
+) {
+    mapboxMap.getStyle { style ->
+        val features = mutableListOf<Feature>()
+
+        for (r in 0 until originalGrid.rows) {
+            for (c in 0 until originalGrid.cols) {
+                val originalElev = originalGrid.getValue(r, c) ?: continue
+
+                // 计算填挖高度
+                val heightDiff = levelingResult.targetElevation - originalElev
+
+                // 获取网格边界
+                val (lon0, lat0, lon1, lat1) = getCellBounds(originalGrid, r, c)
+
+                // 创建多边形要素
+                val ring = listOf(
+                    Point.fromLngLat(lon0, lat0),
+                    Point.fromLngLat(lon1, lat0),
+                    Point.fromLngLat(lon1, lat1),
+                    Point.fromLngLat(lon0, lat1),
+                    Point.fromLngLat(lon0, lat0)
+                )
+                val poly = Polygon.fromLngLats(listOf(ring))
+                val feature = Feature.fromGeometry(poly)
+
+                // 设置属性：根据填挖状态设置颜色
+                when {
+                    heightDiff > 0 -> { // 填方区域
+                        feature.addStringProperty("color", "#FF6B6B") // 红色
+                        feature.addStringProperty("type", "fill")
+                        feature.addNumberProperty("height", heightDiff)
+                    }
+
+                    heightDiff < 0 -> { // 挖方区域
+                        feature.addStringProperty("color", "#4ECDC4") // 青色
+                        feature.addStringProperty("type", "cut")
+                        feature.addNumberProperty("height", abs(heightDiff))
+                    }
+
+                    else -> { // 无变化区域
+                        feature.addStringProperty("color", "#F7FFF7") // 浅绿色
+                        feature.addStringProperty("type", "no-change")
+                        feature.addNumberProperty("height", 0.0)
+                    }
+                }
+
+                features.add(feature)
+            }
+        }
+
+        // 添加图例说明
+        addLegendFeature(features, originalGrid, levelingResult)
+
+        // 设置图层
+        setupEarthworkLayer(style, sourceId, layerId, features)
+    }
+}
+
+/**
+ * 修正的绘制函数 - 确保南北方向正确
+ */
+fun MapView.displayLevelingResultCorrected(
+    originalGrid: GridModel,
+    levelingResult: LevelingResult,
+    sourceId: String = "leveling-result",
+    layerId: String = "leveling-layer",
+    palette: (Double?) -> String
+) {
+    // 先验证方向
+    // debugGridOrientation(originalGrid)
+
+    mapboxMap.getStyle { style ->
+        val features = mutableListOf<Feature>()
+        val maxX = lonToMercX(originalGrid.maxLon)
+        val minX = maxX
+        val maxY = latToMercY(originalGrid.maxLat)
+        val cellMeters = originalGrid.cellSizeMeters
+
+        for (r in 0 until originalGrid.rows) {
+            for (c in 0 until originalGrid.cols) {
+                val originalElev = originalGrid.getValue(r, c) ?: continue
+
+                // 计算填挖高度
+                val heightDiff = levelingResult.targetElevation - originalElev
+
+                // 计算栅格边界
+                val x0 = minX + c * cellMeters
+                val y0 = maxY - r * cellMeters
+                val x1 = x0 + cellMeters
+                val y1 = y0 - cellMeters
+
+                val lon0 = mercXToLon(x0)
+                val lat0 = mercYToLat(y0)
+                val lon1 = mercXToLon(x1)
+                val lat1 = mercYToLat(y1)
+
+                // 创建多边形要素 - 确保顶点顺序正确（逆时针）
+                val ring = listOf(
+                    Point.fromLngLat(lon0, lat0),
+                    Point.fromLngLat(lon1, lat0),
+                    Point.fromLngLat(lon1, lat1),
+                    Point.fromLngLat(lon0, lat1),
+                    Point.fromLngLat(lon0, lat0)
+                )
+
+                val poly = Polygon.fromLngLats(listOf(ring))
+                val feature = Feature.fromGeometry(poly)
+
+                feature.addStringProperty("color", palette(levelingResult.targetElevation))
+
+                if (false) {
+                    // 设置属性
+                    when {
+                        heightDiff > 1.0 -> {
+                            feature.addStringProperty("color", "#FF0000")
+                            feature.addStringProperty("type", "heavy-fill")
+                        }
+
+                        heightDiff > 0 -> {
+                            feature.addStringProperty("color", "#FF6B6B")
+                            feature.addStringProperty("type", "light-fill")
+                        }
+
+                        heightDiff < -1.0 -> {
+                            feature.addStringProperty("color", "#0066CC")
+                            feature.addStringProperty("type", "heavy-cut")
+                        }
+
+                        heightDiff < 0 -> {
+                            feature.addStringProperty("color", "#4ECDC4")
+                            feature.addStringProperty("type", "light-cut")
+                        }
+
+                        else -> {
+                            feature.addStringProperty("color", "#F7FFF7")
+                            feature.addStringProperty("type", "no-change")
+                        }
+                    }
+
+                    feature.addNumberProperty("height_diff", heightDiff)
+                    feature.addNumberProperty("row", r.toDouble())
+                    feature.addNumberProperty("col", c.toDouble())
+                }
+
+                features.add(feature)
+            }
+        }
+
+        // 设置图层
+        setupEarthworkLayer(style, sourceId, layerId, features)
+
+        println("✅ 绘制完成：共 ${features.size} 个网格单元")
+    }
+}
+
+/**
+ * 绘制斜坡设计结果
+ */
+fun MapView.displaySlopeResult(
+    originalGrid: GridModel,
+    slopeResult: SlopeResult,
+    sourceId: String = "slope-result",
+    layerId: String = "slope-layer",
+    palette: (Double?) -> String,
+    showElevationText: Boolean = false
+) {
+    val elevationList = mutableListOf<Double>()
+    mapboxMap.getStyle { style ->
+        val features = mutableListOf<Feature>()
+        val designGrid = slopeResult.designSurface
+
+        // val minX = lonToMercX(originalGrid.minLon)
+        // 对比测试，将绘制到原来图形的左边
+        val minX = lonToMercX(originalGrid.minLon * 2 - originalGrid.maxLon)
+        val maxY = latToMercY(originalGrid.maxLat)
+        val cellMeters = originalGrid.cellSizeMeters
+
+        for (r in 0 until originalGrid.rows) {
+            for (c in 0 until originalGrid.cols) {
+                val originalElev = originalGrid.getValue(r, c) ?: continue
+                val designElev = designGrid.getValue(r, c) ?: continue
+                elevationList.add(designElev)
+
+                // 计算填挖高度
+                val heightDiff = designElev - originalElev
+
+                // 计算栅格边界
+                val x0 = minX + c * cellMeters
+                val y0 = maxY - r * cellMeters
+                val x1 = x0 + cellMeters
+                val y1 = y0 - cellMeters
+
+                val lon0 = mercXToLon(x0)
+                val lat0 = mercYToLat(y0)
+                val lon1 = mercXToLon(x1)
+                val lat1 = mercYToLat(y1)
+
+                // 1. 创建多边形要素（背景色）
+                val ring = listOf(
+                    Point.fromLngLat(lon0, lat0),
+                    Point.fromLngLat(lon1, lat0),
+                    Point.fromLngLat(lon1, lat1),
+                    Point.fromLngLat(lon0, lat1),
+                    Point.fromLngLat(lon0, lat0)
+                )
+                val poly = Polygon.fromLngLats(listOf(ring))
+                val feature = Feature.fromGeometry(poly)
+
+                // 显示高差
+                // feature.addStringProperty("color", palette(heightDiff))
+                // 显示设计高度，测试坡向是否正确，和高度是否计算正确
+                feature.addStringProperty("color", palette(designElev))
+                // 显示原始高度
+                // feature.addStringProperty("color", palette(originalElev))
+                features.add(feature)
+
+                if (showElevationText) {
+                    val centerLon = (lon0 + lon1) / 2
+                    val centerLat = (lat0 + lat1) / 2
+
+                    val textPoint = Point.fromLngLat(centerLon, centerLat)
+                    val textFeature = Feature.fromGeometry(textPoint)
+
+                    textFeature.addStringProperty("text", "%.1f".format(designElev))
+                    textFeature.addStringProperty("type", "elevation-text")
+                    textFeature.addNumberProperty("elevation", designElev)
+                    textFeature.addNumberProperty("row", r.toDouble())
+                    textFeature.addNumberProperty("col", c.toDouble())
+
+                    features.add(textFeature)
+                }
+            }
+        }
+
+        Log.d("displayGridWithDirectionArrows", "对比区域的土方量计算: ${elevationList.sum()}, 平均值：${elevationList.average()}")
+
+        // 添加坡向箭头
+        addSlopeDirectionArrow(features, originalGrid, slopeResult.slopeDirection)
+
+        // 设置图层
+        setupEarthworkLayer(style, sourceId, layerId, features, showElevationText = showElevationText)
+    }
+}
+
+/**
+ * 添加坡向箭头
+ */
+private fun addSlopeDirectionArrow(
+    features: MutableList<Feature>,
+    grid: GridModel,
+    slopeDirection: Double
+) {
+    val centerLon = (grid.minLon + grid.maxLon) / 2
+    val centerLat = (grid.minLat + grid.maxLat) / 2
+
+    // 箭头长度（基于区域大小）
+    val arrowLength = (grid.maxLon - grid.minLon) * 0.2
+
+    val arrowDirectionRad = Math.toRadians(slopeDirection)
+    val endLon = centerLon + cos(arrowDirectionRad) * arrowLength
+    val endLat = centerLat + sin(arrowDirectionRad) * arrowLength
+
+    // 创建箭头线
+    val arrowLine = LineString.fromLngLats(
+        listOf(
+            Point.fromLngLat(centerLon, centerLat),
+            Point.fromLngLat(endLon, endLat)
+        )
+    )
+
+    val arrowFeature = Feature.fromGeometry(arrowLine)
+    arrowFeature.addStringProperty("color", "#000000") // 黑色箭头
+    arrowFeature.addStringProperty("type", "direction-arrow")
+    features.add(arrowFeature)
+}
+
+/**
+ * 添加图例说明
+ */
+private fun addLegendFeature(
+    features: MutableList<Feature>,
+    grid: GridModel,
+    levelingResult: LevelingResult
+) {
+    // 在图左上角添加文本说明
+    val legendLon = grid.minLon + (grid.maxLon - grid.minLon) * 0.02
+    val legendLat = grid.maxLat - (grid.maxLat - grid.minLat) * 0.02
+
+    // 这里可以添加图例要素，但Mapbox对文本支持有限
+    // 可以考虑用点要素+属性，然后在客户端显示弹窗
+    val legendPoint = Point.fromLngLat(legendLon, legendLat)
+    val legendFeature = Feature.fromGeometry(legendPoint)
+    legendFeature.addStringProperty("type", "legend")
+    legendFeature.addStringProperty("title", "填平结果")
+    legendFeature.addStringProperty(
+        "content",
+        "目标高程: ${"%.2f".format(levelingResult.targetElevation)}m\n" +
+                "挖方: ${"%.0f".format(levelingResult.earthworkResult.cutVolume)}m³\n" +
+                "填方: ${"%.0f".format(levelingResult.earthworkResult.fillVolume)}m³"
+    )
+    features.add(legendFeature)
+}
+
+/**
+ * 获取网格单元边界
+ */
+private fun getCellBounds(grid: GridModel, row: Int, col: Int): Quadruple<Double, Double, Double, Double> {
+    val lon0 = grid.minLon + col * (grid.maxLon - grid.minLon) / grid.cols
+    val lon1 = grid.minLon + (col + 1) * (grid.maxLon - grid.minLon) / grid.cols
+    val lat0 = grid.minLat + row * (grid.maxLat - grid.minLat) / grid.rows
+    val lat1 = grid.minLat + (row + 1) * (grid.maxLat - grid.minLat) / grid.rows
+
+    return Quadruple(lon0, lat0, lon1, lat1)
+}
+
+/**
+ * 完整的土方工程图层设置 - 修正版
+ */
+private fun MapView.setupEarthworkLayer(
+    style: Style,
+    sourceId: String,
+    layerId: String,
+    features: List<Feature>,
+    showElevationText: Boolean = false
+) {
+    // 创建数据源
+    val source = geoJsonSource(sourceId) {
+        featureCollection(FeatureCollection.fromFeatures(features))
+    }
+
+    // 清理旧图层
+    try {
+        style.removeStyleLayer(layerId)
+    } catch (_: Exception) {
+    }
+    try {
+        style.removeStyleLayer("$layerId-arrow")
+    } catch (_: Exception) {
+    }
+    try {
+        style.removeStyleLayer("$layerId-outline")
+    } catch (_: Exception) {
+    }
+    try {
+        style.removeStyleLayer("$layerId-text")
+    } catch (_: Exception) {
+    }
+
+    if (style.styleSourceExists(sourceId)) {
+        style.removeStyleSource(sourceId)
+    }
+
+    // 添加数据源
+    style.addSource(source)
+
+    // 主填充图层
+    val fillLayer = FillLayer(layerId, sourceId).apply {
+        fillColor(Expression.toColor(Expression.get("color")))
+        fillOpacity(0.7)
+    }
+    style.addLayer(fillLayer)
+
+    // 边框图层
+    val outlineLayer = LineLayer("$layerId-outline", sourceId).apply {
+        lineColor("#333333")
+        lineWidth(1.0)
+        lineOpacity(0.5)
+    }
+    style.addLayer(outlineLayer)
+
+    // 箭头图层
+    val arrowLayer = LineLayer("$layerId-arrow", sourceId).apply {
+        lineColor(Expression.toColor(Expression.get("color")))
+        lineWidth(3.0)
+        lineCap(LineCap.ROUND)
+        lineJoin(LineJoin.ROUND)
+        filter(
+            Expression.eq(
+                Expression.get("type"),
+                Expression.literal("direction-arrow")
+            )
+        )
+    }
+    style.addLayer(arrowLayer)
+
+    if (showElevationText) {
+        val textLayer = SymbolLayer("$layerId-text", sourceId).apply {
+            textField(Expression.get("text"))
+            textSize(10.0)
+            textColor("#000000")
+            textHaloColor("#FFFFFF")
+            textHaloWidth(1.0)
+            textOpacity(0.9)
+            textAllowOverlap(false)
+            textIgnorePlacement(false)
+            textAnchor(Expression.literal("center"))
+            filter(Expression.eq(Expression.get("type"), Expression.literal("elevation-text")))
+        }
+        style.addLayer(textLayer)
+    }
+}
+
+// 辅助数据类
+data class Quadruple<A, B, C, D>(val first: A, val second: B, val third: C, val fourth: D)

app/src/main/java/com/icegps/geotools/FindContours.kt

@@ -0,0 +1,319 @@
+package com.icegps.geotools
+
+import android.util.Log
+import com.icegps.common.helper.GeoHelper
+import com.icegps.geotools.ktx.TAG
+import com.icegps.geotools.ktx.niceStr
+import com.icegps.math.geometry.Vector3D
+import com.mapbox.geojson.Feature
+import com.mapbox.geojson.FeatureCollection
+import com.mapbox.geojson.LineString
+import com.mapbox.geojson.Point
+import com.mapbox.maps.MapView
+import com.mapbox.maps.Style
+import com.mapbox.maps.extension.style.expressions.generated.Expression
+import com.mapbox.maps.extension.style.layers.addLayer
+import com.mapbox.maps.extension.style.layers.generated.LineLayer
+import com.mapbox.maps.extension.style.sources.addSource
+import com.mapbox.maps.extension.style.sources.generated.geoJsonSource
+import org.openrndr.math.Vector2
+import org.openrndr.shape.LineSegment
+import org.openrndr.shape.ShapeContour
+import kotlin.math.max
+import kotlin.math.min
+
+/**
+ * Find contours for a function [f] using the marching squares algorithm. A contour is found when f(x) crosses zero.
+ * @param useInterpolation intersection points will be interpolated if true, default true
+ * @return a list of [ShapeContour] instances
+ */
+fun findContours(
+    grid: GridModel,
+    useInterpolation: Boolean = true
+): List<ShapeContour> {
+    val segments = mutableListOf<LineSegment>()
+    val segmentsMap = mutableMapOf<Vector2, MutableList<LineSegment>>()
+
+    val corner = Vector2(grid.minLon, grid.maxLat)
+    val targetElevation = 27.0
+    println("cells = ${grid.cells.joinToString { it.toString() }}")
+
+    val zero = 0.0
+    for (y in 0 until grid.rows) {
+        for (x in 0 until grid.cols) {
+
+            // Here we check if we are at a right or top border. This is to ensure we create closed contours
+            // later on in the process.
+
+            fun getValue(x: Int, y: Int): Double? {
+                return (grid.getValue(x, y) ?: return null) - targetElevation
+            }
+
+            val v00 = if (x == 0 || y == 0) zero else (getValue(x, y) ?: zero)
+            val v10 = if (y == 0) zero else (getValue(x + 1, y) ?: zero)
+            val v01 = if (x == 0) zero else (getValue(x, y + 1) ?: zero)
+            val v11 = getValue(x + 1, y + 1) ?: zero
+
+            val p00 = Vector2(x.toDouble(), y.toDouble()) * grid.cellSizeMeters + corner
+            val p10 = Vector2((x + 1).toDouble(), y.toDouble()) * grid.cellSizeMeters + corner
+            val p01 = Vector2(x.toDouble(), (y + 1).toDouble()) * grid.cellSizeMeters + corner
+            val p11 = Vector2((x + 1).toDouble(), (y + 1).toDouble()) * grid.cellSizeMeters + corner
+
+            val v = buildString {
+                append("v00 = ${v00}, ")
+                append("v10 = ${v10}, ")
+                append("v01 = ${v01}, ")
+                append("v11 = ${v11}, ")
+            }
+            println("v = ${v}")
+            val p = buildString {
+                append("p00 = ${p00}, ")
+                append("p10 = ${p10}, ")
+                append("p01 = ${p01}, ")
+                append("p11 = ${p11}, ")
+            }
+            println("p = ${p}")
+
+
+            val index = (if (v00 >= 0.0) 1 else 0) +
+                    (if (v10 >= 0.0) 2 else 0) +
+                    (if (v01 >= 0.0) 4 else 0) +
+                    (if (v11 >= 0.0) 8 else 0)
+
+            println("index = ${index}")
+
+            fun blend(v1: Double, v2: Double): Double {
+                if (useInterpolation) {
+                    require(v1 == v1 && v2 == v2)
+                    val f1 = min(v1, v2)
+                    val f2 = max(v1, v2)
+                    val v = (-f1) / (f2 - f1)
+
+                    require(v == v)
+                    require(v in 0.0..1.0)
+
+                    return if (f1 == v1) {
+                        v
+                    } else {
+                        1.0 - v
+                    }
+                } else {
+                    return 0.5
+                }
+            }
+
+            fun emitLine(
+                p00: Vector2, p01: Vector2, v00: Double, v01: Double,
+                p10: Vector2, p11: Vector2, v10: Double, v11: Double
+            ) {
+                val r0 = blend(v00, v01)
+                val r1 = blend(v10, v11)
+
+                val v0 = p00.mix(p01, r0)
+                val v1 = p10.mix(p11, r1)
+                val l0 = LineSegment(v0, v1)
+                segmentsMap.getOrPut(v1) { mutableListOf() }.add(l0)
+                segmentsMap.getOrPut(v0) { mutableListOf() }.add(l0)
+                segments.add(l0)
+            }
+
+            when (index) {
+                0, 15 -> {}
+                1, 15 xor 1 -> {
+                    emitLine(p00, p01, v00, v01, p00, p10, v00, v10)
+                }
+
+                2, 15 xor 2 -> {
+                    emitLine(p00, p10, v00, v10, p10, p11, v10, v11)
+                }
+
+                3, 15 xor 3 -> {
+                    emitLine(p00, p01, v00, v01, p10, p11, v10, v11)
+                }
+
+                4, 15 xor 4 -> {
+                    emitLine(p00, p01, v00, v01, p01, p11, v01, v11)
+                }
+
+                5, 15 xor 5 -> {
+                    emitLine(p00, p10, v00, v10, p01, p11, v01, v11)
+                }
+
+                6, 15 xor 6 -> {
+                    emitLine(p00, p01, v00, v01, p00, p10, v00, v10)
+                    emitLine(p01, p11, v01, v11, p10, p11, v10, v11)
+                }
+
+                7, 15 xor 7 -> {
+                    emitLine(p01, p11, v01, v11, p10, p11, v10, v11)
+                }
+            }
+        }
+    }
+
+    val processedSegments = mutableSetOf<LineSegment>()
+    val contours = mutableListOf<ShapeContour>()
+    for (segment in segments) {
+        if (segment in processedSegments) {
+            continue
+        } else {
+            val collected = mutableListOf<Vector2>()
+            var current: LineSegment? = segment
+            var closed = true
+            var lastVertex = Vector2.INFINITY
+            do {
+                current!!
+                if (lastVertex.squaredDistanceTo(current.start) > 1E-5) {
+                    collected.add(current.start)
+                }
+                lastVertex = current.start
+                processedSegments.add(current)
+                if (segmentsMap[current.start]!!.size < 2) {
+                    closed = false
+                }
+                val hold = current
+                current = segmentsMap[current.start]?.firstOrNull { it !in processedSegments }
+                if (current == null) {
+                    current = segmentsMap[hold.end]?.firstOrNull { it !in processedSegments }
+                }
+            } while (current != segment && current != null)
+
+            contours.add(ShapeContour.fromPoints(collected, closed = closed))
+        }
+    }
+    return contours
+}
+
+fun MapView.displayFindContours(
+    grid: GridModel,
+    sourceId: String = "findContours-source-id",
+    layerId: String = "findContours-layer-id",
+    palette: (Double?) -> String
+) {
+    val contours = findContours(grid = grid)
+    val features = mutableListOf<Feature>()
+    contours.forEachIndexed { index, shapeContour: ShapeContour ->
+        val points = shapeContour.segments.map {
+            val start = it.start
+            Point.fromLngLat(start.x, start.y)
+        }
+        val lineString = LineString.fromLngLats(points)
+        val feature = Feature.fromGeometry(lineString)
+        val elevation = (index / contours.size) * 255.0
+        feature.addStringProperty("color", "#FF0000")
+        feature.addNumberProperty("elevation", elevation)
+        feature.addStringProperty("type", "contour-line")
+        features.add(feature)
+    }
+    mapboxMap.getStyle { style ->
+        setupContourLayers(
+            style = style,
+            sourceId = sourceId,
+            layerId = layerId,
+            features = features
+        )
+    }
+}
+
+fun MapView.displayFindContours(
+    grid: GridModel,
+    sourceId: String = "findContours-source-id",
+    layerId: String = "findContours-layer-id",
+) {
+    val contours = findContours(grid = grid)
+    mapboxMap.getStyle { style ->
+        val features = mutableListOf<Feature>()
+
+        // 添加调试信息
+        Log.d("ContoursDebug", "找到 ${contours.size} 条等高线")
+
+        contours.forEachIndexed { index, shapeContour ->
+            // 调试每条等高线
+            Log.d("ContoursDebug", "等高线 $index: ${shapeContour.segments.size} 个线段")
+
+            val points = shapeContour.segments.map { segment ->
+                Point.fromLngLat(segment.start.x, segment.start.y)
+            }.distinct()
+            if (points.isNotEmpty())
+            Log.d(
+                TAG,
+                buildString {
+                    appendLine("displayFindContours: ")
+                    val geoHelper = GeoHelper.getSharedInstance()
+                    val wsg84List = points.map {
+                        Vector3D(x = it.longitude(), y = it.latitude(), z = 20.0)
+                    }
+                    appendLine(wsg84List.niceStr())
+                    val enus = points.map {
+                        geoHelper.wgs84ToENU(lon = it.longitude(), lat = it.latitude(), hgt = 20.0)
+                    }.map {
+                        Vector3D(it.x, it.y, it.z)
+                    }
+                    appendLine(enus.niceStr())
+                    appendLine("minX = ${enus.minOf { it.x }}")
+                    appendLine("maxX = ${enus.maxOf { it.x }}")
+                    appendLine("minY = ${enus.minOf { it.y }}")
+                    appendLine("maxY = ${enus.maxOf { it.y }}")
+                }
+            )
+            // 检查坐标范围
+            if (points.isNotEmpty()) {
+                val lngs = points.map { it.longitude() }
+                val lats = points.map { it.latitude() }
+                Log.d("ContoursDebug", "坐标范围: 经度 ${lngs.minOrNull()}~${lngs.maxOrNull()}, 纬度 ${lats.minOrNull()}~${lats.maxOrNull()}")
+            }
+
+            val lineString = LineString.fromLngLats(points)
+            val feature = Feature.fromGeometry(lineString)
+
+            feature.addStringProperty("color", "#FF0000")
+            feature.addNumberProperty("elevation", 1.0)
+            feature.addStringProperty("type", "contour-line")
+
+            features.add(feature)
+        }
+
+        Log.d("ContoursDebug", "总共创建了 ${features.size} 个要素")
+        setupFindContoursLayers(style, sourceId, layerId, features)
+    }
+}
+
+/**
+ * 设置等高线图层
+ */
+fun MapView.setupFindContoursLayers(
+    style: Style,
+    sourceId: String,
+    layerId: String,
+    features: List<Feature>
+) {
+    val source = geoJsonSource(sourceId) {
+        featureCollection(FeatureCollection.fromFeatures(features))
+    }
+
+    // 清理旧图层
+    try {
+        style.removeStyleLayer(layerId)
+    } catch (_: Exception) {
+    }
+    try {
+        style.removeStyleLayer("$layerId-labels")
+    } catch (_: Exception) {
+    }
+
+    if (style.styleSourceExists(sourceId)) {
+        style.removeStyleSource(sourceId)
+    }
+
+    // 添加数据源
+    style.addSource(source)
+
+    // 等高线图层
+    val lineLayer = LineLayer(layerId, sourceId).apply {
+        lineColor(Expression.toColor(Expression.get("color")))
+        lineWidth(2.0)
+        lineOpacity(0.8)
+        filter(Expression.eq(Expression.get("type"), Expression.literal("contour-line")))
+    }
+    style.addLayer(lineLayer)
+}

app/src/main/java/com/icegps/geotools/GeoJsonUtils.kt

@@ -13,7 +13,7 @@ import com.mapbox.geojson.Polygon
  */
 object GeoJsonUtils {
     // 生成三角形（Polygon）FeatureCollection
-    fun <T : IPoint> trianglesToPolygons(delaunator: Delaunator<T>): FeatureCollection {
+    fun <T : IPoint> trianglesToPolygons(delaunator: IDelaunator<T>): FeatureCollection {
         val features = mutableListOf<Feature>()
         val tris = delaunator.triangles
         // triangles 是按 3 个索引为一组三角形存储
@@ -48,7 +48,7 @@ object GeoJsonUtils {
     }
 
     // 生成边（LineString）FeatureCollection（每条边只输出一次）
-    fun <T : IPoint> trianglesToUniqueEdges(delaunator: Delaunator<T>): FeatureCollection {
+    fun <T : IPoint> trianglesToUniqueEdges(delaunator: IDelaunator<T>): FeatureCollection {
         val features = mutableListOf<Feature>()
         val seen = HashSet<Pair<Int, Int>>()
         val tris = delaunator.triangles

app/src/main/java/com/icegps/geotools/GridCell.kt

@@ -10,17 +10,23 @@ import android.graphics.Canvas
 import android.graphics.Color
 import android.graphics.Paint
 import android.graphics.RectF
-import androidx.core.graphics.toColorInt
-import com.icegps.geotools.model.DPoint
+import android.util.Log
+import com.icegps.geotools.model.IGeoPoint
+import com.icegps.math.geometry.Vector2D
 import com.mapbox.geojson.Feature
 import com.mapbox.geojson.FeatureCollection
+import com.mapbox.geojson.LineString
 import com.mapbox.geojson.Point
 import com.mapbox.geojson.Polygon
 import com.mapbox.maps.MapView
 import com.mapbox.maps.extension.style.expressions.generated.Expression
 import com.mapbox.maps.extension.style.layers.addLayer
 import com.mapbox.maps.extension.style.layers.generated.FillLayer
+import com.mapbox.maps.extension.style.layers.generated.LineLayer
+import com.mapbox.maps.extension.style.layers.generated.SymbolLayer
 import com.mapbox.maps.extension.style.layers.generated.rasterLayer
+import com.mapbox.maps.extension.style.layers.properties.generated.LineCap
+import com.mapbox.maps.extension.style.layers.properties.generated.LineJoin
 import com.mapbox.maps.extension.style.layers.properties.generated.Visibility
 import com.mapbox.maps.extension.style.sources.addSource
 import com.mapbox.maps.extension.style.sources.generated.ImageSource
@@ -29,13 +35,17 @@ import com.mapbox.maps.extension.style.sources.generated.imageSource
 import com.mapbox.maps.extension.style.sources.getSourceAs
 import com.mapbox.maps.extension.style.sources.updateImage
 import kotlin.math.PI
+import kotlin.math.abs
 import kotlin.math.absoluteValue
 import kotlin.math.atan
+import kotlin.math.atan2
 import kotlin.math.ceil
+import kotlin.math.cos
 import kotlin.math.exp
 import kotlin.math.ln
 import kotlin.math.max
 import kotlin.math.min
+import kotlin.math.sin
 import kotlin.math.tan
 
 // -----------------------------
@@ -56,10 +66,8 @@ fun mercYToLat(y: Double): Double {
 // -----------------------------
 // Geometry helpers
 // -----------------------------
-data class Vec2(val x: Double, val y: Double)
-
 /** 点是否在三角形内（在 mercator 坐标系中） — 使用重心 / 矩阵法 */
-fun pointInTriangle(pt: Vec2, a: Vec2, b: Vec2, c: Vec2): Boolean {
+fun pointInTriangle(pt: Vector2D, a: Vector2D, b: Vector2D, c: Vector2D): Boolean {
     val v0x = c.x - a.x
     val v0y = c.y - a.y
     val v1x = b.x - a.x
@@ -84,9 +92,9 @@ fun pointInTriangle(pt: Vec2, a: Vec2, b: Vec2, c: Vec2): Boolean {
 /** 可选：用三角形顶点值做双线性/重心内插（这里示例：按顶点值插值）
  * valueAtVerts: DoubleArray of length 3 for the triangle's vertex values
  */
-fun barycentricInterpolate(pt: Vec2, a: Vec2, b: Vec2, c: Vec2, values: DoubleArray): Double {
+fun barycentricInterpolateLegacy(pt: Vector2D, a: Vector2D, b: Vector2D, c: Vector2D, values: DoubleArray): Double {
     // compute areas (using cross product) as barycentric weights
-    val area = { p1: Vec2, p2: Vec2, p3: Vec2 ->
+    val area = { p1: Vector2D, p2: Vector2D, p3: Vector2D ->
         ((p2.x - p1.x) * (p3.y - p1.y) - (p3.x - p1.x) * (p2.y - p1.y)).absoluteValue / 2.0
     }
     val areaTotal = area(a, b, c)
@@ -97,22 +105,95 @@ fun barycentricInterpolate(pt: Vec2, a: Vec2, b: Vec2, c: Vec2, values: DoubleAr
     return values[0] * wA + values[1] * wB + values[2] * wC
 }
 
+fun barycentricInterpolate(pt: Vector2D, a: Vector2D, b: Vector2D, c: Vector2D, values: DoubleArray): Double {
+    val denom = (b.y - c.y) * (a.x - c.x) + (c.x - b.x) * (a.y - c.y)
+
+    // 避免除零（退化三角形）
+    if (abs(denom) < 1e-10) {
+        return values.average()  // 或返回第一个值
+    }
+
+    val w1 = ((b.y - c.y) * (pt.x - c.x) + (c.x - b.x) * (pt.y - c.y)) / denom
+    val w2 = ((c.y - a.y) * (pt.x - c.x) + (a.x - c.x) * (pt.y - c.y)) / denom
+    val w3 = 1.0 - w1 - w2
+
+    // 验证点是否在三角形内（可选）
+    if (w1 < 0 || w2 < 0 || w3 < 0 || w1 > 1 || w2 > 1 || w3 > 1) {
+        // 点在外部的处理策略
+        return when {
+            w1 < -0.1 || w2 < -0.1 || w3 < -0.1 -> Double.NaN  // 太远，返回无效值
+            else -> values[0] * w1 + values[1] * w2 + values[2] * w3  // 轻微外部仍插值
+        }
+    }
+
+    return values[0] * w1 + values[1] * w2 + values[2] * w3
+}
+
 // -----------------------------
 // 主函数：把 Delaunay 转成规则栅格（格子中心采样）
 // -----------------------------
 data class GridCell(val row: Int, val col: Int, val centerLon: Double, val centerLat: Double, var value: Double? = null)
 
+/**
+ * 栅格数据模型
+ * 用于表示规则网格化的地理空间数据，如数字高程模型(DEM)、遥感影像等
+ *
+ * @property minLon 最小经度 (度)，定义栅格区域的西部边界
+ * @property minLat 最小纬度 (度)，定义栅格区域的南部边界
+ * @property maxLon 最大经度 (度)，定义栅格区域的东部边界
+ * @property maxLat 最大纬度 (度)，定义栅格区域的北部边界
+ * @property rows 栅格行数 (Y方向)，表示南北方向的网格数量
+ * @property cols 栅格列数 (X方向)，表示东西方向的网格数量
+ * @property cellSizeMeters 像元大小 (米)，每个网格单元的实际物理尺寸
+ * @property cells 栅格数据数组，存储每个像元的值，如高程、温度、像素值等
+ *                 数组长度: rows * cols
+ *                 存储顺序: 行优先 (row-major)，索引计算: idx = r * cols + c
+ *                 其中 r 为行索引(0 ~ rows-1)，c 为列索引(0 ~ cols-1)
+ *                 允许空值(Double?)表示无数据区域
+ *
+ * 示例:
+ * 一个3x3的栅格，数据布局如下:
+ *
+ *     经度方向 (X/Columns) →
+ * 纬 [ (0,0)  (0,1)  (0,2) ] ← 第0行
+ * 度 [ (1,0)  (1,1)  (1,2) ] ← 第1行
+ * 方 [ (2,0)  (2,1)  (2,2) ] ← 第2行
+ * 向 ↑
+ *
+ * 在cells数组中的存储顺序:
+ * [ (0,0), (0,1), (0,2), (1,0), (1,1), (1,2), (2,0), (2,1), (2,2) ]
+ */
 data class GridModel(
-    val minLon: Double, val minLat: Double,
-    val maxLon: Double, val maxLat: Double,
+    val minLon: Double,
+    val minLat: Double,
+    val maxLon: Double,
+    val maxLat: Double,
     val rows: Int,
     val cols: Int,
     val cellSizeMeters: Double,
     val cells: Array<Double?> // length rows*cols, row-major: idx = r*cols + c
-)
+) {
+    /**
+     * 获取指定行列位置的像元值
+     * @param row 行索引 (0 ~ rows-1)
+     * @param col 列索引 (0 ~ cols-1)
+     * @return 像元值，如果位置无效或为无数据返回null
+     */
+    fun getValue(row: Int, col: Int): Double? {
+        if (row < 0 || row >= rows || col < 0 || col >= cols) {
+            return null
+        }
+        return cells[row * cols + col]
+    }
+
+}
+
+/**
+ * 现在的需求就是，根据某个区域的 GridModel或者三角网，计算当前这块区域的土方量，第一个需求就是填平，用当前计算的土方量计算出平均高度，把这块区域填平，第二个需求就是，给出一个坡向，用这块区域的土做坡
+ */
 
 fun triangulationToGrid(
-    delaunator: Delaunator<DPoint>,
+    delaunator: IDelaunator<IGeoPoint>,
     cellSizeMeters: Double = 50.0, // 每个格子的边长（米）
     maxSidePixels: Int = 5000 // 限制 max rows/cols 防止 OOM（可选）
 ): GridModel {
@@ -158,7 +239,7 @@ fun triangulationToGrid(
     val cells = Array<Double?>(rows * cols) { null }
 
     // 准备点/三角形在 mercator 下的缓存坐标
-    val mercPts = pts.map { p -> Vec2(lonToMercX(p.x), latToMercY(p.y)) }
+    val mercPts = pts.map { p -> Vector2D(lonToMercX(p.x), latToMercY(p.y)) }
 
     // triangles 数组（每 3 个为一组）
     val triIdx = delaunator.triangles
@@ -166,9 +247,10 @@ fun triangulationToGrid(
 
     // For potential vertex values: if you have scalar per vertex, prepare here.
     // Example: create placeholder values (e.g., 0.0). Replace with your actual values if available.
-    val vertexValues = DoubleArray(pts.size) { 0.0 }
 
-    // 3) iterate triangles and rasterize onto grid by checking the grid cells that intersect triangle bbox
+    val vertexValues = pts.map { it.z }
+
+    // 3) 通过检查与三角形 bbox 相交的网格单元来迭代三角形并栅格化到网格上
     for (ti in 0 until triCount) {
         val i0 = triIdx[3 * ti]
         val i1 = triIdx[3 * ti + 1]
@@ -177,19 +259,19 @@ fun triangulationToGrid(
         val b = mercPts[i1]
         val c = mercPts[i2]
 
-        // triangle bbox in mercator
+        // 墨卡托三角形 bbox
         val tminX = minOf(a.x, b.x, c.x)
         val tmaxX = maxOf(a.x, b.x, c.x)
         val tminY = minOf(a.y, b.y, c.y)
         val tmaxY = maxOf(a.y, b.y, c.y)
 
-        // convert bbox to grid indices (clamp)
+        // 将 bbox 转换为网格索引（clamp
         val colMin = ((tminX - minX) / cellSizeMeters).toInt().coerceIn(0, cols - 1)
         val colMax = ((tmaxX - minX) / cellSizeMeters).toInt().coerceIn(0, cols - 1)
         val rowMin = ((maxY - tmaxY) / cellSizeMeters).toInt().coerceIn(0, rows - 1) // 注意 Y 方向
         val rowMax = ((maxY - tminY) / cellSizeMeters).toInt().coerceIn(0, rows - 1)
 
-        // optional: get vertex values for interpolation
+        // 可选：获取插值的顶点值
         val triVertexVals = doubleArrayOf(vertexValues[i0], vertexValues[i1], vertexValues[i2])
 
         for (r in rowMin..rowMax) {
@@ -197,16 +279,17 @@ fun triangulationToGrid(
                 // center of this cell in mercator
                 val centerX = minX + (cIdx + 0.5) * cellSizeMeters
                 val centerY = maxY - (r + 0.5) * cellSizeMeters
-                val pt = Vec2(centerX, centerY)
+                val pt = Vector2D(centerX, centerY)
                 if (pointInTriangle(pt, a, b, c)) {
                     // example: set cell value as triangle index, or do interpolation
                     // cells index:
                     val idx = r * cols + cIdx
                     // choose value: triangle index -> convert to Double
-                    cells[idx] = ti.toDouble()
+                    // cells[idx] = ti.toDouble()
                     // OR for interpolation:
-                    // val valInterp = barycentricInterpolate(pt, a, b, c, triVertexVals)
-                    // cells[idx] = valInterp
+                    val valInterp = barycentricInterpolateLegacy(pt, a, b, c, triVertexVals)
+                    // Log.d("triangulationToGrid", "triangulationToGrid: 从 barycentricInterpolate 计算的 cell[idx] = ${valInterp}")
+                    cells[idx] = valInterp
                 }
             }
         }
@@ -400,6 +483,97 @@ fun MapView.displayGridAsImageSourceHighRes(
     }
 }
 
+class SimplePalette(
+    private var range: ClosedFloatingPointRange<Double>
+) {
+    fun setRange(range: ClosedFloatingPointRange<Double>) {
+        this.range = range
+    }
+
+    private val colors: Map<Int, String>
+
+    init {
+        colors = generateTerrainColorMap()
+    }
+
+    fun palette(value: Double?): String {
+        if (value == null) return "#00000000"
+        val minH = range.start
+        val maxH = range.endInclusive
+        val normalized = ((value - minH) / (maxH - minH)).coerceIn(0.0, 1.0)
+        return colors[(normalized * 255).toInt()] ?: "#00000000"
+    }
+
+    fun palette1(value: Double?): String {
+        return if (value == null) "#00000000" else {
+            // 假设您已经知道高度范围，或者动态计算
+            val minH = range.start
+            val maxH = range.endInclusive
+            val normalized = ((value - minH) / (maxH - minH)).coerceIn(0.0, 1.0)
+            val alpha = (normalized * 255).toInt()
+            String.format("#%02X%02X%02X", alpha, 0, 0)
+        }.also {
+            Log.d("simplePalette", "$value -> $it")
+        }
+    }
+
+    fun generateTerrainColorMap(): MutableMap<Int, String> {
+        val colorMap = mutableMapOf<Int, String>()
+
+        // 定义关键颜色点
+        val blue = Color(0, 0, 255)    // 低地势 - 蓝色
+        val cyan = Color(0, 255, 255)  // 中间过渡
+        val green = Color(0, 255, 0)   // 中间过渡
+        val yellow = Color(255, 255, 0) // 中间过渡
+        val red = Color(255, 0, 0)     // 高地势 - 红色
+
+        for (i in 0..255) {
+            val position = i / 255.0
+
+            val color = when {
+                position < 0.25 -> interpolateColor(blue, cyan, position / 0.25)
+                position < 0.5 -> interpolateColor(cyan, green, (position - 0.25) / 0.25)
+                position < 0.75 -> interpolateColor(green, yellow, (position - 0.5) / 0.25)
+                else -> interpolateColor(yellow, red, (position - 0.75) / 0.25)
+            }
+            colorMap[i] = color.toHex()
+        }
+
+        return colorMap
+    }
+
+    fun interpolateColor(start: Color, end: Color, fraction: Double): Color {
+        val r = (start.red + (end.red - start.red) * fraction).toInt()
+        val g = (start.green + (end.green - start.green) * fraction).toInt()
+        val b = (start.blue + (end.blue - start.blue) * fraction).toInt()
+        return Color(r, g, b)
+    }
+
+    // Color类简化实现
+    class Color(val red: Int, val green: Int, val blue: Int) {
+        fun toArgb(): Int {
+            return (0xFF shl 24) or (red shl 16) or (green shl 8) or blue
+        }
+
+        fun toHex(): String {
+            return String.format("#%06X", toArgb() and 0xFFFFFF)
+        }
+    }
+}
+
+// 简化的使用方式（如果您想要最简洁的）
+val simplePalette: (Double?) -> String = { value ->
+    if (value == null) "#00000000" else {
+        // 假设您已经知道高度范围，或者动态计算
+        val minH = 0.0
+        val maxH = 255.0
+        val normalized = ((value - minH) / (maxH - minH)).coerceIn(0.0, 1.0)
+        val alpha = (normalized * 255).toInt()
+        String.format("#%02X%02X%02X%02X", alpha, alpha, alpha, alpha)
+    }.also {
+        Log.d("simplePalette", "${value} -> ${it}")
+    }
+}
 
 // -----------------------------
 // 可选：把每个格子做成 GeoJSON Polygons（每格一个 Fill）并显示（交互式，但格子多时非常慢）
@@ -466,29 +640,109 @@ fun MapView.displayGridAsGeoJsonPolygons(
     }
 }
 
-fun MapView.displayGridAsGeoJsonWithHeight(
+/**
+ * 计算栅格单元的地势降低方向（下坡方向）
+ * 返回角度（0-360度，0表示北，90表示东）
+ */
+fun GridModel.calculateSlopeDirection(row: Int, col: Int): Double? {
+    // 获取当前单元和相邻单元的高程
+    val center = getValue(row, col) ?: return null
+
+    // 注意：在墨卡托坐标系中，Y轴正向是北，负向是南
+    // 但在我们的栅格索引中，row=0 是最北边，row增大向南
+    val north = getValueSafe(row - 1, col)  // 上一行（北边）
+    val south = getValueSafe(row + 1, col)  // 下一行（南边）
+    val east = getValueSafe(row, col + 1)   // 右边（东边）
+    val west = getValueSafe(row, col - 1)   // 左边（西边）
+
+    // 需要至少两个有效相邻点才能计算方向
+    val validNeighbors = listOfNotNull(north, south, east, west)
+    if (validNeighbors.size < 2) return null
+
+    // 计算梯度 - 修正Y轴方向
+    val dzdx = when {
+        east != null && west != null -> (east - west) / (2 * cellSizeMeters)
+        east != null -> (east - center) / cellSizeMeters
+        west != null -> (center - west) / cellSizeMeters
+        else -> 0.0
+    }
+
+    // 关键修正：Y轴方向
+    // 在墨卡托中，Y增大是向北，但在我们的栅格索引中，row增大是向南
+    // 所以dy应该是 (north - south)，因为north在更小的row索引上
+    val dzdy = when {
+        north != null && south != null -> (north - south) / (2 * cellSizeMeters)  // 修正：北减南
+        north != null -> (north - center) / cellSizeMeters
+        south != null -> (center - south) / cellSizeMeters  // 注意符号
+        else -> 0.0
+    }
+
+    println("调试: row=$row, col=$col, center=$center, north=$north, south=$south, east=$east, west=$west")
+    println("调试: dzdx=$dzdx, dzdy=$dzdy")
+
+    // 计算坡向（地势升高的方向）
+    var aspect = atan2(dzdy, dzdx) * 180 / Math.PI
+    if (aspect < 0) aspect += 360.0
+
+    println("调试: 坡向=$aspect")
+
+    // 返回地势降低的方向（坡向的相反方向）
+    val downhill = (aspect + 180.0) % 360.0
+    println("调试: 下坡方向=$downhill")
+
+    return downhill
+}
+
+/**
+ * 安全的获取相邻单元值
+ */
+private fun GridModel.getValueSafe(row: Int, col: Int): Double? {
+    return if (row in 0 until rows && col in 0 until cols) {
+        getValue(row, col)
+    } else {
+        null
+    }
+}
+
+fun MapView.displayGridWithDirectionArrows(
     grid: GridModel,
-    testSourceId: String,
-    testLayerId: String,
-    heightToColor: (Double) -> Int
+    polygonSourceId: String,
+    polygonLayerId: String,
+    arrowSourceId: String,
+    arrowLayerId: String,
+    palette: (Double?) -> String,
+    arrowScale: Double = 0.3,
+    minSlopeThreshold: Double = 0.01,
+    arrowEnabled: Boolean = false,
+    showElevationText: Boolean = false
 ) {
+    val elevationList = mutableListOf<Double>()
     mapboxMap.getStyle { style ->
+        val polygonFeatures = mutableListOf<Feature>()
+        val arrowFeatures = mutableListOf<Feature>()
+        val textFeatures = mutableListOf<Feature>()
 
-        val features = mutableListOf<Feature>()
-
-        // 计算每格经纬度跨度
-        val deltaLon = (grid.maxLon - grid.minLon) / grid.cols
-        val deltaLat = (grid.maxLat - grid.minLat) / grid.rows
+        val minX = lonToMercX(grid.minLon)
+        val maxY = latToMercY(grid.maxLat)
+        val cellMeters = grid.cellSizeMeters
 
         for (r in 0 until grid.rows) {
             for (c in 0 until grid.cols) {
-                val z = grid.cells[r * grid.cols + c] ?: continue
+                val idx = r * grid.cols + c
+                val v = grid.cells[idx] ?: continue
 
-                val lon0 = grid.minLon + c * deltaLon
-                val lon1 = grid.minLon + (c + 1) * deltaLon
-                val lat0 = grid.maxLat - r * deltaLat
-                val lat1 = grid.maxLat - (r + 1) * deltaLat
+                // 计算栅格边界
+                val x0 = minX + c * cellMeters
+                val y0 = maxY - r * cellMeters
+                val x1 = x0 + cellMeters
+                val y1 = y0 - cellMeters
 
+                val lon0 = mercXToLon(x0)
+                val lat0 = mercYToLat(y0)
+                val lon1 = mercXToLon(x1)
+                val lat1 = mercYToLat(y1)
+
+                // 1. 创建多边形要素（背景色）
                 val ring = listOf(
                     Point.fromLngLat(lon0, lat0),
                     Point.fromLngLat(lon1, lat0),
@@ -496,34 +750,172 @@ fun MapView.displayGridAsGeoJsonWithHeight(
                     Point.fromLngLat(lon0, lat1),
                     Point.fromLngLat(lon0, lat0)
                 )
-
                 val poly = Polygon.fromLngLats(listOf(ring))
-                val f = Feature.fromGeometry(poly)
+                val polyFeature = Feature.fromGeometry(poly)
+                polyFeature.addStringProperty("color", palette(v))
+                polyFeature.addNumberProperty("value", v ?: -9999.0)
+                polygonFeatures.add(polyFeature)
 
-                // 添加高度属性
-                f.addNumberProperty("value", z)
+                // 2. 创建箭头要素（如果可计算方向）
+                val slopeDirection = if (arrowEnabled) grid.calculateSlopeDirection(r, c) else null
+                if (slopeDirection != null) {
+                    // 计算箭头位置（栅格中心）
+                    val centerLon = (lon0 + lon1) / 2
+                    val centerLat = (lat0 + lat1) / 2
 
-                // 根据回调生成颜色
-                val colorInt = heightToColor(z)
-                val colorStr = String.format("#%08X", colorInt)
-                f.addStringProperty("color", colorStr)
+                    // 计算箭头长度（栅格尺寸的一部分）
+                    val arrowLength = cellMeters * arrowScale
 
-                features.add(f)
+                    // 计算箭头终点（地势降低的方向）
+                    val arrowDirectionRad = Math.toRadians(slopeDirection)
+                    val endX = x0 + cellMeters / 2 + cos(arrowDirectionRad) * arrowLength
+                    val endY = y0 - cellMeters / 2 + sin(arrowDirectionRad) * arrowLength
+
+                    val endLon = mercXToLon(endX)
+                    val endLat = mercYToLat(endY)
+
+                    // 创建线要素表示箭杆
+                    val arrowLine = LineString.fromLngLats(
+                        listOf(
+                            Point.fromLngLat(centerLon, centerLat),
+                            Point.fromLngLat(endLon, endLat)
+                        )
+                    )
+
+                    val arrowFeature = Feature.fromGeometry(arrowLine)
+                    arrowFeature.addNumberProperty("direction", slopeDirection)
+                    arrowFeature.addStringProperty("color", "#FF0000")
+                    arrowFeatures.add(arrowFeature)
+
+                    // 添加三角形箭头头部
+                    val headSize = cellMeters * 0.15
+                    val leftRad = arrowDirectionRad + Math.PI * 0.8
+                    val rightRad = arrowDirectionRad - Math.PI * 0.8
+
+                    val leftX = endX + cos(leftRad) * headSize
+                    val leftY = endY + sin(leftRad) * headSize
+                    val rightX = endX + cos(rightRad) * headSize
+                    val rightY = endY + sin(rightRad) * headSize
+
+                    val leftLon = mercXToLon(leftX)
+                    val leftLat = mercYToLat(leftY)
+                    val rightLon = mercXToLon(rightX)
+                    val rightLat = mercYToLat(rightY)
+
+                    // 创建三角形箭头头部
+                    val headRing = listOf(
+                        Point.fromLngLat(endLon, endLat),
+                        Point.fromLngLat(leftLon, leftLat),
+                        Point.fromLngLat(rightLon, rightLat),
+                        Point.fromLngLat(endLon, endLat)
+                    )
+                    val headPolygon = Polygon.fromLngLats(listOf(headRing))
+                    val headFeature = Feature.fromGeometry(headPolygon)
+                    headFeature.addStringProperty("color", "#FF0000")
+                    arrowFeatures.add(headFeature)
+                }
+                if (showElevationText) {
+                    val originalElev = grid.getValue(r, c) ?: continue
+                    elevationList.add(originalElev)
+                    val centerLon = (lon0 + lon1) / 2
+                    val centerLat = (lat0 + lat1) / 2
+
+                    val textPoint = Point.fromLngLat(centerLon, centerLat)
+                    val textFeature = Feature.fromGeometry(textPoint)
+
+                    textFeature.addStringProperty("text", "%.1f".format(originalElev))
+                    textFeature.addStringProperty("type", "elevation-text")
+                    textFeature.addNumberProperty("elevation", originalElev)
+                    textFeature.addNumberProperty("row", r.toDouble())
+                    textFeature.addNumberProperty("col", c.toDouble())
+
+                    textFeatures.add(textFeature)
+                }
             }
         }
 
-        val fc = FeatureCollection.fromFeatures(features)
+        Log.d("displayGridWithDirectionArrows", "运来区域的土方量计算: ${elevationList.sum()}, 平均值：${elevationList.average()}")
 
-        // 添加或更新 GeoJSON Source
-        if (style.styleSourceExists(testSourceId)) style.removeStyleSource(testSourceId)
-        style.addSource(geoJsonSource(testSourceId) { featureCollection(fc) })
+        // 3. 处理多边形图层 - 先移除Layer再移除Source
+        try {
+            style.removeStyleLayer(polygonLayerId)
+        } catch (_: Exception) {
+        }
 
-        // 创建 FillLayer 并使用 feature.color
-        try { style.removeStyleLayer(testLayerId) } catch (_: Exception) {}
-        val fillLayer = FillLayer(testLayerId, testSourceId).apply {
+        if (style.styleSourceExists(polygonSourceId)) {
+            style.removeStyleSource(polygonSourceId)
+        }
+
+        val polygonSource = geoJsonSource(polygonSourceId) {
+            featureCollection(FeatureCollection.fromFeatures(polygonFeatures))
+        }
+        style.addSource(polygonSource)
+
+        val fillLayer = FillLayer(polygonLayerId, polygonSourceId).apply {
             fillColor(Expression.toColor(Expression.get("color")))
-            fillOpacity(0.9)
+            fillOpacity(0.7)
         }
         style.addLayer(fillLayer)
+
+        // 4. 处理箭头图层 - 先移除Layer再移除Source
+        val arrowHeadLayerId = "$arrowLayerId-head"
+
+        try {
+            style.removeStyleLayer(arrowLayerId)
+        } catch (_: Exception) {
+        }
+
+        try {
+            style.removeStyleLayer(arrowHeadLayerId)
+        } catch (_: Exception) {
+        }
+
+        if (style.styleSourceExists(arrowSourceId)) {
+            style.removeStyleSource(arrowSourceId)
+        }
+
+        val arrowSource = geoJsonSource(arrowSourceId) {
+            featureCollection(FeatureCollection.fromFeatures(arrowFeatures))
+        }
+        style.addSource(arrowSource)
+
+        val lineLayer = LineLayer(arrowLayerId, arrowSourceId).apply {
+            lineColor(Expression.toColor(Expression.get("color")))
+            lineWidth(2.0)
+            lineCap(LineCap.ROUND)
+            lineJoin(LineJoin.ROUND)
+        }
+        style.addLayer(lineLayer)
+
+        val headLayer = FillLayer(arrowHeadLayerId, arrowSourceId).apply {
+            fillColor(Expression.toColor(Expression.get("color")))
+            fillOpacity(1.0)
+        }
+        style.addLayer(headLayer)
+
+        // 绘制数字
+        val textSourceId = "text-source-id-0"
+        val textLayerId = "text-layer-id-0"
+        style.removeStyleLayer(textLayerId)
+        if (style.styleSourceExists(textSourceId)) {
+            style.removeStyleSource(textSourceId)
+        }
+        val textSource = geoJsonSource(textSourceId) {
+            featureCollection(FeatureCollection.fromFeatures(textFeatures))
+        }
+        style.addSource(textSource)
+        val textLayer = SymbolLayer(textLayerId, textSourceId).apply {
+            textField(Expression.get("text"))
+            textSize(10.0)
+            textColor("#000000")
+            textHaloColor("#FFFFFF")
+            textHaloWidth(1.0)
+            textOpacity(0.9)
+            textAllowOverlap(false)
+            textIgnorePlacement(false)
+            textAnchor(Expression.literal("center"))
+            filter(Expression.eq(Expression.get("type"), Expression.literal("elevation-text")))
+        }
+        style.addLayer(textLayer)
     }
 }

app/src/main/java/com/icegps/geotools/HeightmapVolcanoGenerator.kt

@@ -0,0 +1,375 @@
+package com.icegps.geotools
+
+/**
+ * @author tabidachinokaze
+ * @date 2025/11/19
+ */
+import com.icegps.math.geometry.Vector3D
+import kotlin.math.PI
+import kotlin.math.abs
+import kotlin.math.cos
+import kotlin.math.exp
+import kotlin.math.floor
+import kotlin.math.max
+import kotlin.math.sin
+import kotlin.math.sqrt
+
+object HeightmapVolcanoGenerator {
+
+    // 基础火山高度图
+    fun generateVolcanoHeightmap(
+        width: Int = 100,
+        height: Int = 100,
+        centerX: Double = 50.0,
+        centerY: Double = 50.0,
+        maxHeight: Double = 60.0,
+        craterRadius: Double = 8.0,
+        volcanoRadius: Double = 30.0
+    ): List<Vector3D> {
+        val points = mutableListOf<Vector3D>()
+
+        for (x in 0 until width) {
+            for (y in 0 until height) {
+                // 计算到火山中心的距离
+                val dx = x - centerX
+                val dy = y - centerY
+                val distance = sqrt(dx * dx + dy * dy)
+
+                // 计算基础火山高度
+                var z = calculateVolcanoHeight(distance, craterRadius, volcanoRadius, maxHeight)
+
+                // 添加噪声细节
+                val noise = perlinNoise(x * 0.1, y * 0.1, 0.1) * 3.0
+                z = max(0.0, z + noise)
+
+                points.add(Vector3D(x.toDouble(), y.toDouble(), z))
+            }
+        }
+
+        return points
+    }
+
+    // 复合火山群高度图
+    fun generateVolcanoClusterHeightmap(
+        width: Int = 150,
+        height: Int = 150,
+        volcanoCount: Int = 3
+    ): List<Vector3D> {
+        val points = mutableListOf<Vector3D>()
+        val volcanoes = generateRandomVolcanoPositions(volcanoCount, width, height)
+
+        for (x in 0 until width) {
+            for (y in 0 until height) {
+                var totalZ = 0.0
+
+                // 叠加所有火山的影响
+                for (volcano in volcanoes) {
+                    val dx = x - volcano.x
+                    val dy = y - volcano.y
+                    val distance = sqrt(dx * dx + dy * dy)
+
+                    if (distance <= volcano.radius) {
+                        val volcanoHeight = calculateVolcanoHeight(
+                            distance,
+                            volcano.craterRadius,
+                            volcano.radius,
+                            volcano.maxHeight
+                        )
+                        totalZ += volcanoHeight
+                    }
+                }
+
+                // 基础地形
+                val baseNoise = perlinNoise(x * 0.02, y * 0.02, 0.05) * 5.0
+                val detailNoise = perlinNoise(x * 0.1, y * 0.1, 0.2) * 2.0
+
+                points.add(Vector3D(x.toDouble(), y.toDouble(), totalZ + baseNoise + detailNoise))
+            }
+        }
+
+        return points
+    }
+
+    // 带熔岩流的火山高度图
+    fun generateVolcanoWithLavaHeightmap(
+        width: Int = 100,
+        height: Int = 100
+    ): List<Vector3D> {
+        val points = mutableListOf<Vector3D>()
+        val centerX = width / 2.0
+        val centerY = height / 2.0
+
+        // 生成熔岩流路径
+        val lavaFlows = generateLavaFlowPaths(centerX, centerY, 3)
+
+        for (x in 0 until width) {
+            for (y in 0 until height) {
+                val dx = x - centerX
+                val dy = y - centerY
+                val distance = sqrt(dx * dx + dy * dy)
+
+                // 基础火山高度
+                var z = calculateVolcanoHeight(distance, 10.0, 35.0, 70.0)
+
+                // 添加熔岩流
+                z += calculateLavaFlowEffect(x.toDouble(), y.toDouble(), lavaFlows)
+
+                // 侵蚀效果
+                z += calculateErosionEffect(x.toDouble(), y.toDouble(), distance, z)
+
+                points.add(Vector3D(x.toDouble(), y.toDouble(), max(0.0, z)))
+            }
+        }
+
+        return points
+    }
+
+    // 破火山口高度图
+    fun generateCalderaHeightmap(
+        width: Int = 100,
+        height: Int = 100
+    ): List<Vector3D> {
+        val points = mutableListOf<Vector3D>()
+        val centerX = width / 2.0
+        val centerY = height / 2.0
+
+        for (x in 0 until width) {
+            for (y in 0 until height) {
+                val dx = x - centerX
+                val dy = y - centerY
+                val distance = sqrt(dx * dx + dy * dy)
+
+                var z = calculateCalderaHeight(distance, 15.0, 45.0, 50.0)
+
+                // 内部平坦区域细节
+                if (distance < 20) {
+                    z += perlinNoise(x * 0.2, y * 0.2, 0.3) * 1.5
+                }
+
+                points.add(Vector3D(x.toDouble(), y.toDouble(), max(0.0, z)))
+            }
+        }
+
+        return points
+    }
+
+    // 线性火山链高度图
+    fun generateVolcanoChainHeightmap(
+        width: Int = 200,
+        height: Int = 100
+    ): List<Vector3D> {
+        val points = mutableListOf<Vector3D>()
+
+        // 在一条线上生成多个火山
+        val chainCenters = listOf(
+            Vector3D(30.0, 50.0, 0.0),
+            Vector3D(70.0, 50.0, 0.0),
+            Vector3D(110.0, 50.0, 0.0),
+            Vector3D(150.0, 50.0, 0.0),
+            Vector3D(170.0, 50.0, 0.0)
+        )
+
+        for (x in 0 until width) {
+            for (y in 0 until height) {
+                var totalZ = 0.0
+
+                for (center in chainCenters) {
+                    val dx = x - center.x
+                    val dy = y - center.y
+                    val distance = sqrt(dx * dx + dy * dy)
+
+                    if (distance <= 25.0) {
+                        val volcanoZ = calculateVolcanoHeight(distance, 6.0, 25.0, 40.0)
+                        totalZ += volcanoZ
+                    }
+                }
+
+                // 添加基底地形，模拟山脉链
+                val baseRidge = calculateMountainRidge(x.toDouble(), y.toDouble(), width, height)
+                totalZ += baseRidge
+
+                points.add(Vector3D(x.toDouble(), y.toDouble(), totalZ))
+            }
+        }
+
+        return points
+    }
+
+    // 辅助函数
+    private data class VolcanoInfo(
+        val x: Double,
+        val y: Double,
+        val radius: Double,
+        val craterRadius: Double,
+        val maxHeight: Double
+    )
+
+    private data class LavaFlowInfo(
+        val startX: Double,
+        val startY: Double,
+        val angle: Double, // 弧度
+        val length: Double,
+        val width: Double,
+        val intensity: Double
+    )
+
+    private fun calculateVolcanoHeight(
+        distance: Double,
+        craterRadius: Double,
+        volcanoRadius: Double,
+        maxHeight: Double
+    ): Double {
+        return when {
+            distance <= craterRadius -> {
+                // 火山口 - 中心凹陷
+                val craterDepth = maxHeight * 0.4
+                craterDepth * (1.0 - distance / craterRadius)
+            }
+
+            distance <= volcanoRadius -> {
+                // 火山锥
+                val slopeDistance = distance - craterRadius
+                val maxSlopeDistance = volcanoRadius - craterRadius
+                val normalized = slopeDistance / maxSlopeDistance
+                maxHeight * (1.0 - normalized * normalized)
+            }
+
+            else -> 0.0
+        }
+    }
+
+    private fun calculateCalderaHeight(
+        distance: Double,
+        innerRadius: Double,
+        outerRadius: Double,
+        rimHeight: Double
+    ): Double {
+        return when {
+            distance <= innerRadius -> {
+                // 平坦的破火山口底部
+                rimHeight * 0.2
+            }
+
+            distance <= outerRadius -> {
+                // 陡峭的边缘
+                val rimDistance = distance - innerRadius
+                val rimWidth = outerRadius - innerRadius
+                val normalized = rimDistance / rimWidth
+                rimHeight * (1.0 - (1.0 - normalized) * (1.0 - normalized))
+            }
+
+            else -> {
+                // 外部平缓斜坡
+                val externalDistance = distance - outerRadius
+                rimHeight * exp(-externalDistance * 0.08)
+            }
+        }
+    }
+
+    private fun calculateLavaFlowEffect(x: Double, y: Double, lavaFlows: List<LavaFlowInfo>): Double {
+        var effect = 0.0
+
+        for (flow in lavaFlows) {
+            val dx = x - flow.startX
+            val dy = y - flow.startY
+
+            // 计算到熔岩流中心线的距离
+            val flowDirX = cos(flow.angle)
+            val flowDirY = sin(flow.angle)
+
+            val projection = dx * flowDirX + dy * flowDirY
+
+            if (projection in 0.0..flow.length) {
+                val perpendicularX = dx - projection * flowDirX
+                val perpendicularY = dy - projection * flowDirY
+                val perpendicularDist = sqrt(perpendicularX * perpendicularX + perpendicularY * perpendicularY)
+
+                if (perpendicularDist <= flow.width) {
+                    val widthFactor = 1.0 - (perpendicularDist / flow.width)
+                    val lengthFactor = 1.0 - (projection / flow.length)
+                    effect += flow.intensity * widthFactor * lengthFactor
+                }
+            }
+        }
+
+        return effect
+    }
+
+    private fun calculateErosionEffect(x: Double, y: Double, distance: Double, height: Double): Double {
+        // 基于坡度的侵蚀
+        val slopeNoise = perlinNoise(x * 0.15, y * 0.15, 0.1) * 2.0
+        // 基于距离的侵蚀
+        val distanceErosion = if (distance > 25) perlinNoise(x * 0.08, y * 0.08, 0.05) * 1.5 else 0.0
+        return slopeNoise + distanceErosion
+    }
+
+    private fun calculateMountainRidge(x: Double, y: Double, width: Int, height: Int): Double {
+        // 创建山脉基底
+        val ridgeCenter = height / 2.0
+        val distanceToRidge = abs(y - ridgeCenter)
+        val ridgeWidth = height * 0.3
+
+        if (distanceToRidge <= ridgeWidth) {
+            val ridgeFactor = 1.0 - (distanceToRidge / ridgeWidth)
+            return ridgeFactor * 15.0 * perlinNoise(x * 0.01, y * 0.01, 0.02)
+        }
+        return 0.0
+    }
+
+    private fun generateRandomVolcanoPositions(count: Int, width: Int, height: Int): List<VolcanoInfo> {
+        return List(count) {
+            VolcanoInfo(
+                x = (width * 0.2 + random() * width * 0.6),
+                y = (height * 0.2 + random() * height * 0.6),
+                radius = 20.0 + random() * 20.0,
+                craterRadius = 5.0 + random() * 7.0,
+                maxHeight = 25.0 + random() * 35.0
+            )
+        }
+    }
+
+    private fun generateLavaFlowPaths(centerX: Double, centerY: Double, count: Int): List<LavaFlowInfo> {
+        return List(count) {
+            LavaFlowInfo(
+                startX = centerX,
+                startY = centerY,
+                angle = random() * 2 * PI,
+                length = 20.0 + random() * 15.0,
+                width = 2.0 + random() * 3.0,
+                intensity = 5.0 + random() * 8.0
+            )
+        }
+    }
+
+    private fun perlinNoise(x: Double, y: Double, frequency: Double): Double {
+        // 简化的柏林噪声实现
+        val x0 = floor(x * frequency)
+        val y0 = floor(y * frequency)
+        val x1 = x0 + 1
+        val y1 = y0 + 1
+
+        fun grad(ix: Int, iy: Int): Double {
+            val random = sin(ix * 12.9898 + iy * 78.233) * 43758.5453
+            return (random % 1.0) * 2 - 1
+        }
+
+        fun interpolate(a: Double, b: Double, w: Double): Double {
+            return a + (b - a) * (w * w * (3 - 2 * w))
+        }
+
+        val g00 = grad(x0.toInt(), y0.toInt())
+        val g10 = grad(x1.toInt(), y0.toInt())
+        val g01 = grad(x0.toInt(), y1.toInt())
+        val g11 = grad(x1.toInt(), y1.toInt())
+
+        val tx = x * frequency - x0
+        val ty = y * frequency - y0
+
+        val n0 = interpolate(g00, g10, tx)
+        val n1 = interpolate(g01, g11, tx)
+
+        return interpolate(n0, n1, ty)
+    }
+
+    private fun random(): Double = Math.random()
+}

app/src/main/java/com/icegps/geotools/MainActivity.kt

@@ -8,16 +8,29 @@ import android.graphics.Canvas
 import android.graphics.Color
 import android.graphics.Paint
 import android.graphics.Path
+import android.graphics.PointF
 import android.graphics.drawable.BitmapDrawable
 import android.os.Bundle
 import android.util.Log
-import androidx.appcompat.app.AppCompatActivity
+import androidx.activity.ComponentActivity
 import androidx.core.content.ContextCompat
+import androidx.lifecycle.lifecycleScope
+import com.google.android.material.slider.Slider
 import com.icegps.common.helper.GeoHelper
+import com.icegps.geotools.api.OpenElevation
+import com.icegps.geotools.api.OpenElevationApi
+import com.icegps.geotools.databinding.ActivityMainBinding
 import com.icegps.geotools.ktx.TAG
 import com.icegps.geotools.ktx.niceStr
+import com.icegps.geotools.ktx.toast
 import com.icegps.geotools.model.DPoint
+import com.icegps.geotools.model.GeoPoint
+import com.icegps.geotools.model.IGeoPoint
+import com.icegps.math.geometry.Angle
+import com.icegps.math.geometry.Vector2D
 import com.icegps.math.geometry.Vector3D
+import com.icegps.math.geometry.degrees
+import com.mapbox.android.gestures.MoveGestureDetector
 import com.mapbox.geojson.Feature
 import com.mapbox.geojson.FeatureCollection
 import com.mapbox.geojson.MultiPoint
@@ -25,6 +38,7 @@ import com.mapbox.geojson.Point
 import com.mapbox.geojson.Polygon
 import com.mapbox.maps.CameraOptions
 import com.mapbox.maps.MapView
+import com.mapbox.maps.ScreenCoordinate
 import com.mapbox.maps.Style
 import com.mapbox.maps.extension.style.expressions.generated.Expression.Companion.get
 import com.mapbox.maps.extension.style.layers.addLayer
@@ -47,21 +61,86 @@ import com.mapbox.maps.extension.style.sources.generated.imageSource
 import com.mapbox.maps.extension.style.sources.getSourceAs
 import com.mapbox.maps.extension.style.sources.updateImage
 import com.mapbox.maps.extension.style.style
+import com.mapbox.maps.plugin.gestures.OnMoveListener
+import com.mapbox.maps.plugin.gestures.addOnMapClickListener
+import com.mapbox.maps.plugin.gestures.addOnMoveListener
+import kotlinx.coroutines.CoroutineExceptionHandler
+import kotlinx.coroutines.flow.MutableStateFlow
+import kotlinx.coroutines.flow.combine
+import kotlinx.coroutines.flow.filter
+import kotlinx.coroutines.flow.filterNotNull
+import kotlinx.coroutines.flow.launchIn
+import kotlinx.coroutines.flow.onEach
+import kotlinx.coroutines.flow.update
+import kotlinx.coroutines.launch
+import kotlin.math.cos
+import kotlin.math.log2
 import kotlin.math.roundToInt
+import kotlin.math.sin
 import kotlin.random.Random
 
 
-class MainActivity : AppCompatActivity() {
+class MainActivity : ComponentActivity() {
     private lateinit var mapView: MapView
     private val geoHelper = GeoHelper.getSharedInstance()
+    private lateinit var binding: ActivityMainBinding
+    private var cellSizeMeters: Double = 100.0
+    private lateinit var delaunator: IDelaunator<IGeoPoint>
+    private val palette = SimplePalette(0.0..255.0)
+    private val delaunatorFlow = MutableStateFlow<IDelaunator<IGeoPoint>?>(null)
+    private val clickedPoints = MutableStateFlow<List<IGeoPoint>>(emptyList())
+    private val openElevation: OpenElevationApi = OpenElevation(
+        client = SharedHttpClient(json = SharedJson())
+    )
+    private val earthworkManager: EarthworkManager = EarthworkManager()
+    private var slopeDirection = MutableStateFlow(0.degrees)
+    private var slopePercentage = MutableStateFlow(50f)
+    private var baseHeightOffset = MutableStateFlow(0f)
+    private var gridModel: GridModel? = null
+    private var centerArrow = MutableStateFlow(Point.fromLngLat(0.0, 0.0))
+    private var arrow = MutableStateFlow(Point.fromLngLat(0.0, 1.0))
+    private var headArrow = MutableStateFlow(emptyList<Point>())
 
     override fun onCreate(savedInstanceState: Bundle?) {
         super.onCreate(savedInstanceState)
-        mapView = MapView(this)
+        //mapView = MapView(this)
 
-        setContentView(mapView)
+        binding = ActivityMainBinding.inflate(layoutInflater, null, false)
+        mapView = binding.mapView
+        setContentView(binding.root)
+        binding.slider.addOnSliderTouchListener(
+            object : Slider.OnSliderTouchListener {
+                override fun onStartTrackingTouch(slider: Slider) {
+                }
+
+                override fun onStopTrackingTouch(slider: Slider) {
+                    val gridModel = triangulationToGrid(
+                        delaunator = delaunator,
+                        cellSizeMeters = cellSizeMeters,
+                        maxSidePixels = 6553500
+                    )
+                    displayDirectionArrows(gridModel)
+                }
+            }
+        )
+        binding.slider.addOnChangeListener { slider, value, fromUser ->
+            cellSizeMeters = value / 100.0
+        }
         initGeoHelper()
-        val pointsVector = coordinateGenerate()
+        // val pointsVector = coordinateGenerate1()
+        val generator = HeightmapVolcanoGenerator
+        val singleVolcano = generator.generateVolcanoHeightmap(
+            width = 100, height = 100,
+            centerX = 50.0, centerY = 50.0,
+            maxHeight = 60.0, craterRadius = 8.0, volcanoRadius = 30.0
+        )
+
+        val volcanoCluster = generator.generateVolcanoClusterHeightmap(100, 100, 4)
+        val lavaVolcano = generator.generateVolcanoWithLavaHeightmap(100, 100)
+        val caldera = generator.generateCalderaHeightmap(100, 100)
+        val volcanoChain = generator.generateVolcanoChainHeightmap(200, 100)
+        // 俯视
+        val pointsVector = volcanoCluster
         val pointsLatLng = pointsVector.map {
             val wgS84 = geoHelper.enuToWGS84Object(
                 GeoHelper.ENU(x = it.x, y = it.y, z = it.z)
@@ -81,11 +160,15 @@ class MainActivity : AppCompatActivity() {
         }*/
         // mapView.mapboxMap.addLayer(symbolLayer)
 
-        val delaunator = Delaunator(points = pointsLatLng.map {
+        delaunator = Delaunator(points = pointsLatLng.map {
             DPoint(x = it.longitude(), y = it.latitude(), z = it.altitude())
         })
         val polygons = GeoJsonUtils.trianglesToPolygons(delaunator)
+        val minZ = delaunator.points.minOf { it.z }
+        val maxZ = delaunator.points.maxOf { it.z }
+        palette.setRange(minZ..maxZ)
         // 显示三角网
+        // 这里可以不需要了，通过 Flow 来加载数据
         if (false) mapView.loadGeoJson(delaunator)
         // val rasterLayer = RasterLayer("raster-layer", sourceId = "raster-source")
         if (false) mapView.loadRaterFromTin(delaunator)
@@ -142,22 +225,28 @@ class MainActivity : AppCompatActivity() {
             cellSizeMeters = 2.0,
             maxSidePixels = 6553500
         )
-        val palette: (Double?) -> Int = { v ->
+
+        /*val palette: (Double?) -> Int = { v ->
             if (v == null) Color.MAGENTA else { // null 显 magenta 便于确认
                 val idx = v.toInt()
                 Color.rgb((idx * 37) and 0xFF, (idx * 61) and 0xFF, (idx * 97) and 0xFF)
             }
-        }
+        }*/
+
         if (false) mapView.displayGridAsImageSourceHighRes(
             grid = gridModel,
             testSourceId = "raster-source-id-1",
             testLayerId = "raster-layer-id-1",
         )
-        if (true) mapView.displayGridAsGeoJsonPolygons(
+
+        // 显示栅格
+        if (false) mapView.displayGridAsGeoJsonPolygons(
             grid = gridModel,
             testSourceId = "raster-source-id-0",
             testLayerId = "raster-layer-id-0",
+            palette = palette::palette
         )
+        displayDirectionArrows(gridModel)
 
         mapView.mapboxMap.setCamera(
             CameraOptions.Builder()
@@ -169,6 +258,260 @@ class MainActivity : AppCompatActivity() {
                 .bearing(0.0)
                 .build()
         )
+        mapView.mapboxMap.addOnMapClickListener { point ->
+            Log.d(TAG, "onCreate: ")
+            lifecycleScope.launch(CoroutineExceptionHandler { coroutineContext, throwable ->
+                toast(throwable.message ?: "请求高程失败")
+            }) {
+                val response = openElevation.lookup(listOf(GeoPoint(longitude = point.longitude(), latitude = point.latitude(), elevation = point.altitude())))
+                val element = response.results.firstOrNull() ?: return@launch
+                toast("[${element.x},${element.y},${element.z}]")
+                clickedPoints.update {
+                    it.toMutableList().apply {
+                        add(element)
+                    }
+                }
+            }
+            true
+        }
+        binding.buttonClear.setOnClickListener {
+            clickedPoints.value = emptyList()
+        }
+        binding.angleSlider.addOnSliderTouchListener(
+            object : Slider.OnSliderTouchListener {
+                override fun onStartTrackingTouch(slider: Slider) {
+                }
+
+                override fun onStopTrackingTouch(slider: Slider) {
+                    this@MainActivity.slopeDirection.value = slider.value.degrees
+                }
+            }
+        )
+        binding.slopePercentageSlider.addOnSliderTouchListener(
+            object : Slider.OnSliderTouchListener {
+                override fun onStartTrackingTouch(slider: Slider) {
+                }
+
+                override fun onStopTrackingTouch(slider: Slider) {
+                    slopePercentage.value = slider.value
+                }
+            }
+        )
+        binding.baseHeightOffsetSlider.addOnSliderTouchListener(
+            object : Slider.OnSliderTouchListener {
+                override fun onStartTrackingTouch(slider: Slider) {
+                }
+
+                override fun onStopTrackingTouch(slider: Slider) {
+                    baseHeightOffset.value = slider.value
+                }
+            }
+        )
+        mapView.mapboxMap.addOnMoveListener(
+            object : OnMoveListener {
+                private var beginning: Boolean = false
+                private var isDragging: Boolean = false
+                private fun getCoordinate(focalPoint: PointF): Point {
+                    return mapView.mapboxMap.coordinateForPixel(ScreenCoordinate(focalPoint.x.toDouble(), focalPoint.y.toDouble()))
+                }
+
+                override fun onMove(detector: MoveGestureDetector): Boolean {
+                    val focalPoint = detector.focalPoint
+                    val point = mapView.mapboxMap.coordinateForPixel(ScreenCoordinate(focalPoint.x.toDouble(), focalPoint.y.toDouble()))
+                    val isPointInPolygon = RayCastingAlgorithm.isPointInPolygon(
+                        point = point,
+                        polygon = headArrow.value
+                    )
+                    Log.d(
+                        TAG,
+                        buildString {
+                            appendLine("onMove: ")
+                            appendLine("${with(focalPoint) { "[$x, $y]" }} -> [${point.longitude()},${point.latitude()}]")
+                            appendLine("headArrow = $isPointInPolygon")
+                        }
+                    )
+                    if (isPointInPolygon) {
+                        isDragging = true
+                    }
+                    if (isDragging) {
+                        arrow.value = point
+                    }
+                    return isDragging
+                }
+
+                override fun onMoveBegin(detector: MoveGestureDetector) {
+                    Log.d(TAG, "onMoveBegin: $detector")
+                    beginning = true
+                }
+
+                override fun onMoveEnd(detector: MoveGestureDetector) {
+                    Log.d(TAG, "onMoveEnd: $detector")
+                    val point = getCoordinate(detector.focalPoint)
+                    val arrow = Vector2D(point.longitude(), point.latitude())
+                    val centerArrowValue = centerArrow.value
+                    val center = Vector2D(centerArrowValue.longitude(), centerArrowValue.latitude())
+                    if (beginning && isDragging) {
+                        slopeDirection.value = (arrow - center).angle
+                    }
+
+                    Log.d(
+                        TAG,
+                        buildString {
+                            appendLine("onMoveEnd: ")
+                            appendLine("${point.longitude()}, ${point.latitude()}")
+                        }
+                    )
+                    isDragging = false
+                    beginning = false
+                }
+            }
+        )
+        initData()
+    }
+
+    private fun displayDirectionArrows(gridModel: GridModel) {
+        this@MainActivity.gridModel = gridModel
+        if (true) mapView.displayGridWithDirectionArrows(
+            grid = gridModel,
+            polygonSourceId = "raster-source-id-2",
+            polygonLayerId = "raster-layer-id-2",
+            arrowSourceId = "arrow-source-id-2",
+            arrowLayerId = "arrow-layer-id-2",
+            palette = palette::palette,
+            arrowScale = 0.3,
+        )
+        fun calculateArrow() {
+            val centerLon = (gridModel.minLon + gridModel.maxLon) / 2
+            val centerLat = (gridModel.minLat + gridModel.maxLat) / 2
+            centerArrow.value = Point.fromLngLat(centerLon, centerLat)
+        }
+        calculateArrow()
+        val trend = gridModel.calculateDominantSlopeTrend()
+        mapView.displayControllableArrow(
+            grid = gridModel,
+            angle = trend.direction.degrees,
+            onHeadArrowChange = {
+                headArrow.value = it
+            }
+        )
+        if (false) mapView.displayOverallTrendArrow(
+            grid = gridModel,
+            trendSourceId = "trend-source-id-0",
+            trendLayerId = "trend-layer-id-0",
+            arrowScale = 0.9,
+            onArrowCenterChange = {
+                centerArrow.value = it
+            },
+            onHeadArrowChange = {
+                headArrow.value = it
+            }
+        )
+        if (false) runCatching {
+            val levelingResult = earthworkManager.calculateLeveling(grid = gridModel)
+
+            // 在Mapbox上显示结果
+            mapView.displayLevelingResultCorrected(
+                gridModel,
+                levelingResult,
+                "leveling-result",
+                "leveling-layer",
+                palette::palette
+            )
+        }.onFailure {
+            it.printStackTrace()
+        }
+        displaySlope(
+            gridModel = gridModel,
+            slopeDirection = slopeDirection.value,
+            slopePercentage = slopePercentage.value.toDouble(),
+            baseHeightOffset = baseHeightOffset.value.toDouble()
+        )
+        if (false) mapView.displayFindContours(
+            grid = gridModel,
+            // palette = palette::palette
+        )
+        if (true) mapView.displayContinuousContours(
+            grid = gridModel,
+        )
+        if (false) mapView.displayContoursByMarchingSquares(
+            grid = gridModel,
+            palette = palette::palette
+        )
+    }
+
+    private fun displaySlope(
+        gridModel: GridModel,
+        slopeDirection: Angle,
+        slopePercentage: Double,
+        baseHeightOffset: Double
+    ) {
+        val slopeResult: SlopeResult = earthworkManager.calculateSlope(
+            grid = gridModel,
+            slopeDirection = slopeDirection.degrees,
+            slopePercentage = slopePercentage,
+            baseHeightOffset = baseHeightOffset
+        )
+
+        Log.d(TAG, slopeResult.earthworkResult.toString())
+        // binding.designCut
+        mapView.displaySlopeResult(
+            gridModel,
+            slopeResult,
+            "slope-result",
+            "slope-layer",
+            palette::palette
+        )
+    }
+
+    private fun initData() {
+        clickedPoints.filter {
+            it.size >= 3
+        }.onEach {
+            delaunator = Delaunator(it)
+            delaunatorFlow.value = delaunator
+        }.launchIn(lifecycleScope)
+        delaunatorFlow.filterNotNull().onEach {
+            mapView.loadGeoJson(it)
+            val minZ = it.points.minOf { it.z }
+            val maxZ = it.points.maxOf { it.z }
+            palette.setRange(minZ..maxZ)
+            val gridModel = triangulationToGrid(
+                delaunator = it,
+                cellSizeMeters = cellSizeMeters,
+                maxSidePixels = 10000000
+            )
+            displayDirectionArrows(gridModel)
+        }.launchIn(lifecycleScope)
+        combine(
+            slopeDirection,
+            slopePercentage,
+            baseHeightOffset
+        ) { slopeDirection, slopePercentage, baseHeightOffset ->
+            gridModel?.let { gridModel ->
+                displaySlope(gridModel, slopeDirection, slopePercentage.toDouble(), baseHeightOffset.toDouble())
+            }
+        }.launchIn(lifecycleScope)
+        combine(
+            centerArrow,
+            arrow
+        ) { center, arrow ->
+            val center = Vector2D(center.longitude(), center.latitude())
+            val arrow = Vector2D(arrow.longitude(), arrow.latitude())
+            val direction = (arrow - center)
+            val atan2 = Angle.atan2(direction.x, direction.y, Vector2D.UP)
+            // val angle = (Angle.FULL - atan2).normalized
+            val angle = atan2.normalized
+            binding.angleSlider.value = angle.degrees.toFloat()
+            gridModel?.let { gridModel ->
+                mapView.displayControllableArrow(
+                    grid = gridModel,
+                    angle = angle,
+                    onHeadArrowChange = {
+                        headArrow.value = it
+                    }
+                )
+            }
+        }.launchIn(lifecycleScope)
     }
 }
 
@@ -182,22 +525,81 @@ fun initGeoHelper(base: DPoint = DPoint(114.476060, 22.771073, 30.897)) {
 }
 
 fun coordinateGenerate(): List<Vector3D> {
+    val minX = -20.0
+    val maxX = 20.0
+    val minY = -20.0
+    val maxY = 20.0
+    val minZ = -20.0
+    val maxZ = 20.0
+    val x: () -> Double = { Random.nextDouble(minX, maxX) }
+    val y: () -> Double = { Random.nextDouble(minY, maxY) }
+    val z: () -> Double = { Random.nextDouble(minZ, maxZ) }
     val dPoints = (0..60).map {
-        Vector3D(x, y, z)
+        Vector3D(x(), y(), z())
     }
     return dPoints
 }
 
-val minX = -20.0
-val maxX = 20.0
-val minY = -20.0
-val maxY = 20.0
-val minZ = -20.0
-val maxZ = 20.0
+fun coordinateGenerate1(): List<Vector3D> {
+    /**
+     * 绕 Z 轴旋转指定角度（弧度）
+     */
+    fun Vector3D.rotateAroundZ(angle: Angle): Vector3D {
+        val cosAngle = cos(angle.radians)
+        val sinAngle = sin(angle.radians)
 
-val x: Double get() = Random.nextDouble(minX, maxX)
-val y: Double get() = Random.nextDouble(minY, maxY)
-val z: Double get() = Random.nextDouble(minZ, maxZ)
+        return Vector3D(
+            x = x * cosAngle - y * sinAngle,
+            y = x * sinAngle + y * cosAngle,
+            z = z
+        )
+    }
+
+    val center = Vector3D()
+    val direction = Vector3D(0.0, 1.0, -1.0)
+    return (0..360).step(10).map {
+        val nowDirection = direction.rotateAroundZ(it.degrees)
+        (0..10).map {
+            center + nowDirection * it
+        }
+    }.flatten()
+}
+
+fun coordinateGenerate3(): List<Vector3D> {
+    val lists = (0..10).map { x ->
+        (0..10).map { y ->
+            Vector3D(x, y, x * y)
+        }
+    }.flatten()
+    return lists
+}
+
+fun coordinateGenerate4(): List<Vector3D> {
+    val center = Vector3D()
+    val forwardDirection = Vector3D.FORWARD
+    val rightDirection = Vector3D.RIGHT
+    val doubles = (0..100).map { log2(it.toDouble()) }
+
+    return (0..50).map { x ->
+        (0..50).map { y ->
+            Vector3D(x.toDouble(), y.toDouble(), log2(x * y * 2.0).coerceAtLeast(0.0))
+        }
+    }.flatten()
+}
+
+fun coordinateGenerateSimple(): List<Vector3D> {
+    val forward = Vector3D.FORWARD
+    val right = Vector3D.RIGHT
+    val up = Vector3D.UP
+
+    return listOf(
+        Vector3D(0.0, 0.0, 0.0),
+        forward * 10 + right * 10 + up * 100,
+        forward * 10 - right * 10,
+        -forward * 10 - right * 10,
+        -forward * 10 + right * 10,
+    )
+}
 
 fun Context.getBitmapFromDrawable(drawableResId: Int): Bitmap {
     val drawable = ContextCompat.getDrawable(this, drawableResId)!!
@@ -219,21 +621,27 @@ fun Context.getBitmapFromDrawable(drawableResId: Int): Bitmap {
 }
 
 // 在 Activity/Fragment 中
-fun MapView.loadGeoJson(delaunator: Delaunator<DPoint>) {
+fun MapView.loadGeoJson(delaunator: IDelaunator<IGeoPoint>) {
+    val lineLayerId = "triangles-line-layer"
+    val sourceId = "triangles-source-id"
+    val fillLayerId = "triangles-fill-layer"
+    mapboxMap.removeStyleLayer(lineLayerId)
+    mapboxMap.removeStyleLayer(fillLayerId)
+    mapboxMap.removeStyleSource(sourceId)
     mapboxMap.getStyle { style ->
         // 1) 构建 FeatureCollection（选择边或三角形）
 
         val edgesFc = GeoJsonUtils.trianglesToUniqueEdges(delaunator) // 或 trianglesToPolygons(delaunator)
 
         // 2) 添加 GeoJsonSource（id = "triangles-source")
-        val source = GeoJsonSource.Builder("triangles-source")
+        val source = GeoJsonSource.Builder(sourceId)
             .featureCollection(edgesFc)
             .build()
         style.addSource(source)
 
         // 3) 添加 LineLayer 渲染三角网边
 
-        val line = lineLayer("triangles-line-layer", "triangles-source") {
+        val line = lineLayer(lineLayerId, sourceId) {
             lineWidth(1.5)
             lineJoin(LineJoin.ROUND)
             lineOpacity(1.0)
@@ -243,11 +651,11 @@ fun MapView.loadGeoJson(delaunator: Delaunator<DPoint>) {
         style.addLayer(line)
 
         // 可选：如果你用 polygons 并想填充三角形
-        val fill = fillLayer("triangles-fill-layer", "triangles-source") {
+        val fill = fillLayer(fillLayerId, sourceId) {
             fillOpacity(0.2)
             fillColor("#00FF00")
         }
-        style.addLayerBelow(fill, "triangles-line-layer")
+        style.addLayerBelow(fill, lineLayerId)
     }
 }
 
@@ -275,7 +683,7 @@ fun MapView.loadImageTest(bbox: BoundingBox) {
     }
 }
 
-fun MapView.loadRaterFromTin(delaunator: Delaunator<DPoint>) {
+fun MapView.loadRaterFromTin(delaunator: IDelaunator<IGeoPoint>) {
     // add ImageSource
     delaunator.points
     val bbox = RasterUtils.boundingBox(delaunator.points)
@@ -290,10 +698,10 @@ fun MapView.loadRaterFromTin(delaunator: Delaunator<DPoint>) {
         valueGetter = { it.z }
     )
     val bboxCoordinates: List<List<Double>> = listOf(
-        listOf(minX, maxY), // top-left  (lon, lat)
-        listOf(maxX, maxY), // top-right
-        listOf(maxX, minY), // bottom-right
-        listOf(minX, minY)  // bottom-left
+        // listOf(minX, maxY), // top-left  (lon, lat)
+        // listOf(maxX, maxY), // top-right
+        // listOf(maxX, minY), // bottom-right
+        // listOf(minX, minY)  // bottom-left
     )
 
     mapboxMap.getStyle { style ->
@@ -419,7 +827,7 @@ fun MapView.loadRasterFromResource() {
 
 @Deprecated("显示效果不对")
 fun MapView.showVoronoiAsPolygons(
-    delaunator: Delaunator<DPoint>,
+    delaunator: IDelaunator<IGeoPoint>,
     testSourceId: String,
     testLayerId: String
 ) {
@@ -460,7 +868,7 @@ fun MapView.showVoronoiAsPolygons(
  *  - pixelsPerDegree: 控制生成的栅格分辨率（每经度多少像素），或直接给 width/height
  */
 fun MapView.rasterizeDelaunayToMap(
-    delaunator: Delaunator<DPoint>,
+    delaunator: IDelaunator<IGeoPoint>,
     testSourceId: String,
     testLayerId: String,
     pixelsPerDegree: Double = 400.0 // 可调：值越大图片越精细、越大

app/src/main/java/com/icegps/geotools/MarchingSquares.kt

@@ -0,0 +1,321 @@
+package com.icegps.geotools
+
+/**
+ * @author tabidachinokaze
+ * @date 2025/11/21
+ */
+import android.util.Log
+import com.icegps.geotools.ktx.TAG
+import com.icegps.math.geometry.Vector2D
+import com.mapbox.geojson.Feature
+import com.mapbox.geojson.LineString
+import com.mapbox.geojson.Point
+import com.mapbox.maps.MapView
+
+/**
+ * 等值线顶点
+ * @param x 经度坐标 (度)
+ * @param y 纬度坐标 (度)
+ */
+typealias Vertex = Vector2D
+
+/**
+ * 等值线 - 由一系列顶点组成的多边形
+ */
+typealias Contour = List<Vertex>
+
+/**
+ * Marching Squares 算法实现
+ * 用于从栅格数据中提取等值线
+ */
+object MarchingSquares {
+    /**
+     * 从栅格数据中提取等值线
+     *
+     * @param grid 栅格数据模型
+     * @param isoValue 等值线对应的数值
+     * @return 等值线列表，每个等值线是一个多边形顶点列表
+     */
+    fun extractContours(grid: GridModel, isoValue: Double): List<Contour> {
+        val contours = mutableListOf<Contour>()
+
+        // 遍历所有网格单元 (注意边界: 0..rows-2, 0..cols-2)
+        for (row in 0 until grid.rows - 1) {
+            for (col in 0 until grid.cols - 1) {
+                val contourSegment = processGridCell(grid, isoValue, row, col)
+                if (contourSegment.isNotEmpty()) {
+                    // 这里简化处理：直接将每个单元的线段作为独立的等值线
+                    // 实际应用中需要将相邻单元的线段连接成完整的等值线
+                    contours.add(contourSegment)
+                }
+            }
+        }
+
+        return connectContourSegments(contours)
+    }
+
+    /**
+     * 处理单个网格单元，生成等值线段
+     *
+     * @param grid 栅格数据
+     * @param isoValue 等值
+     * @param row 单元起始行
+     * @param col 单元起始列
+     * @return 等值线段顶点列表 (0, 1, 或 2个顶点对)
+     */
+    private fun processGridCell(
+        grid: GridModel,
+        isoValue: Double,
+        row: Int,
+        col: Int
+    ): Contour {
+        // 获取单元四个角点的值
+        val topLeft = grid.getValue(row, col)
+        val topRight = grid.getValue(row, col + 1)
+        val bottomLeft = grid.getValue(row + 1, col)
+        val bottomRight = grid.getValue(row + 1, col + 1)
+
+        // 检查是否有无数据点，如果有则跳过该单元
+        if (topLeft == null || topRight == null || bottomLeft == null || bottomRight == null) {
+            return emptyList()
+        }
+
+        // 计算状态索引 (使用标准的Marching Squares权重)
+        var state = 0
+        if (topLeft >= isoValue) state = state or 1      // 左下角: 权重1
+        if (topRight >= isoValue) state = state or 2     // 右下角: 权重2
+        if (bottomRight >= isoValue) state = state or 4  // 右上角: 权重4
+        if (bottomLeft >= isoValue) state = state or 8   // 左上角: 权重8
+
+        Log.d(TAG, "processGridCell: state = ${state}")
+
+        // 根据状态索引确定等值线配置
+        return when (state) {
+            0, 15 -> emptyList() // 全内或全外，无等值线
+            1 -> listOf(
+                interpolateVertex(grid, row, col, row + 1, col, topLeft, bottomLeft, isoValue), // 左边
+                interpolateVertex(grid, row, col, row, col + 1, topLeft, topRight, isoValue)    // 上边
+            )
+
+            2 -> listOf(
+                interpolateVertex(grid, row, col, row, col + 1, topLeft, topRight, isoValue),   // 上边
+                interpolateVertex(grid, row, col + 1, row + 1, col + 1, topRight, bottomRight, isoValue) // 右边
+            )
+
+            3 -> listOf(
+                interpolateVertex(grid, row, col, row + 1, col, topLeft, bottomLeft, isoValue), // 左边
+                interpolateVertex(grid, row, col + 1, row + 1, col + 1, topRight, bottomRight, isoValue) // 右边
+            )
+
+            4 -> listOf(
+                interpolateVertex(grid, row, col + 1, row + 1, col + 1, topRight, bottomRight, isoValue), // 右边
+                interpolateVertex(grid, row + 1, col, row + 1, col + 1, bottomLeft, bottomRight, isoValue) // 下边
+            )
+
+            5 -> {
+                // 歧义情况：使用平均值法解决
+                val centerValue = (topLeft + topRight + bottomLeft + bottomRight) / 4.0
+                if (centerValue >= isoValue) {
+                    listOf(
+                        interpolateVertex(grid, row, col, row + 1, col, topLeft, bottomLeft, isoValue), // 左边
+                        interpolateVertex(grid, row, col + 1, row + 1, col + 1, topRight, bottomRight, isoValue), // 右边
+                        interpolateVertex(grid, row, col, row, col + 1, topLeft, topRight, isoValue),   // 上边
+                        interpolateVertex(grid, row + 1, col, row + 1, col + 1, bottomLeft, bottomRight, isoValue) // 下边
+                    )
+                } else {
+                    listOf(
+                        interpolateVertex(grid, row, col, row, col + 1, topLeft, topRight, isoValue),   // 上边
+                        interpolateVertex(grid, row, col, row + 1, col, topLeft, bottomLeft, isoValue), // 左边
+                        interpolateVertex(grid, row, col + 1, row + 1, col + 1, topRight, bottomRight, isoValue), // 右边
+                        interpolateVertex(grid, row + 1, col, row + 1, col + 1, bottomLeft, bottomRight, isoValue) // 下边
+                    )
+                }
+            }
+
+            6 -> listOf(
+                interpolateVertex(grid, row, col, row, col + 1, topLeft, topRight, isoValue),   // 上边
+                interpolateVertex(grid, row + 1, col, row + 1, col + 1, bottomLeft, bottomRight, isoValue) // 下边
+            )
+
+            7 -> listOf(
+                interpolateVertex(grid, row, col, row + 1, col, topLeft, bottomLeft, isoValue), // 左边
+                interpolateVertex(grid, row + 1, col, row + 1, col + 1, bottomLeft, bottomRight, isoValue) // 下边
+            )
+
+            8 -> listOf(
+                interpolateVertex(grid, row + 1, col, row + 1, col + 1, bottomLeft, bottomRight, isoValue), // 下边
+                interpolateVertex(grid, row, col, row + 1, col, topLeft, bottomLeft, isoValue) // 左边
+            )
+
+            9 -> listOf(
+                interpolateVertex(grid, row, col, row, col + 1, topLeft, topRight, isoValue),   // 上边
+                interpolateVertex(grid, row + 1, col, row + 1, col + 1, bottomLeft, bottomRight, isoValue) // 下边
+            )
+
+            10 -> {
+                // 歧义情况：使用平均值法解决
+                val centerValue = (topLeft + topRight + bottomLeft + bottomRight) / 4.0
+                if (centerValue >= isoValue) {
+                    listOf(
+                        interpolateVertex(grid, row, col, row, col + 1, topLeft, topRight, isoValue),   // 上边
+                        interpolateVertex(grid, row, col, row + 1, col, topLeft, bottomLeft, isoValue), // 左边
+                        interpolateVertex(grid, row, col + 1, row + 1, col + 1, topRight, bottomRight, isoValue), // 右边
+                        interpolateVertex(grid, row + 1, col, row + 1, col + 1, bottomLeft, bottomRight, isoValue) // 下边
+                    )
+                } else {
+                    listOf(
+                        interpolateVertex(grid, row, col, row + 1, col, topLeft, bottomLeft, isoValue), // 左边
+                        interpolateVertex(grid, row, col + 1, row + 1, col + 1, topRight, bottomRight, isoValue), // 右边
+                        interpolateVertex(grid, row, col, row, col + 1, topLeft, topRight, isoValue),   // 上边
+                        interpolateVertex(grid, row + 1, col, row + 1, col + 1, bottomLeft, bottomRight, isoValue) // 下边
+                    )
+                }
+            }
+
+            11 -> listOf(
+                interpolateVertex(grid, row, col + 1, row + 1, col + 1, topRight, bottomRight, isoValue), // 右边
+                interpolateVertex(grid, row + 1, col, row + 1, col + 1, bottomLeft, bottomRight, isoValue) // 下边
+            )
+
+            12 -> listOf(
+                interpolateVertex(grid, row, col, row + 1, col, topLeft, bottomLeft, isoValue), // 左边
+                interpolateVertex(grid, row, col, row, col + 1, topLeft, topRight, isoValue)    // 上边
+            )
+
+            13 -> listOf(
+                interpolateVertex(grid, row, col + 1, row + 1, col + 1, topRight, bottomRight, isoValue), // 右边
+                interpolateVertex(grid, row, col, row, col + 1, topLeft, topRight, isoValue)    // 上边
+            )
+
+            14 -> listOf(
+                interpolateVertex(grid, row, col, row + 1, col, topLeft, bottomLeft, isoValue), // 左边
+                interpolateVertex(grid, row, col + 1, row + 1, col + 1, topRight, bottomRight, isoValue) // 右边
+            )
+
+            else -> emptyList()
+        }
+    }
+
+    /**
+     * 线性插值计算等值线顶点坐标
+     *
+     * @param grid 栅格数据
+     * @param row1 第一个点的行
+     * @param col1 第一个点的列
+     * @param row2 第二个点的行
+     * @param col2 第二个点的列
+     * @param value1 第一个点的值
+     * @param value2 第二个点的值
+     * @param isoValue 等值
+     * @return 插值得到的顶点坐标
+     */
+    private fun interpolateVertex(
+        grid: GridModel,
+        row1: Int, col1: Int,
+        row2: Int, col2: Int,
+        value1: Double, value2: Double,
+        isoValue: Double
+    ): Vertex {
+        // 计算插值比例
+        val t = (isoValue - value1) / (value2 - value1)
+
+        // 计算两个端点的地理坐标
+        val x1 = grid.minLon + col1 * (grid.maxLon - grid.minLon) / (grid.cols - 1)
+        val y1 = grid.minLat + row1 * (grid.maxLat - grid.minLat) / (grid.rows - 1)
+        val x2 = grid.minLon + col2 * (grid.maxLon - grid.minLon) / (grid.cols - 1)
+        val y2 = grid.minLat + row2 * (grid.maxLat - grid.minLat) / (grid.rows - 1)
+
+        // 线性插值计算顶点坐标
+        val x = x1 + t * (x2 - x1)
+        val y = y1 + t * (y2 - y1)
+
+        return Vertex(x, y)
+    }
+
+    /**
+     * 连接分散的等值线段形成完整的等值线
+     * 这是一个简化的实现，实际应用可能需要更复杂的连接逻辑
+     */
+    private fun connectContourSegments(segments: List<Contour>): List<Contour> {
+        if (segments.isEmpty()) return emptyList()
+
+        val connectedContours = mutableListOf<Contour>()
+        val used = BooleanArray(segments.size) { false }
+
+        for (i in segments.indices) {
+            if (!used[i] && segments[i].size >= 2) {
+                val currentContour = mutableListOf<Vertex>()
+                currentContour.addAll(segments[i])
+                used[i] = true
+
+                // 简化的连接逻辑：查找可以连接的线段
+                var foundConnection = true
+                while (foundConnection) {
+                    foundConnection = false
+                    for (j in segments.indices) {
+                        if (!used[j] && segments[j].size >= 2) {
+                            // 检查是否可以连接到当前等值线的首尾
+                            if (isClose(segments[j].first(), currentContour.last())) {
+                                currentContour.addAll(segments[j].drop(1))
+                                used[j] = true
+                                foundConnection = true
+                                break
+                            } else if (isClose(segments[j].last(), currentContour.first())) {
+                                currentContour.addAll(0, segments[j].dropLast(1))
+                                used[j] = true
+                                foundConnection = true
+                                break
+                            }
+                        }
+                    }
+                }
+
+                if (currentContour.size >= 2) {
+                    connectedContours.add(currentContour)
+                }
+            }
+        }
+
+        return connectedContours
+    }
+
+    /**
+     * 判断两个顶点是否足够接近（用于线段连接）
+     */
+    private fun isClose(v1: Vertex, v2: Vertex, tolerance: Double = 1e-6): Boolean {
+        return Math.abs(v1.x - v2.x) < tolerance && Math.abs(v1.y - v2.y) < tolerance
+    }
+}
+
+fun MapView.displayContoursByMarchingSquares(
+    grid: GridModel,
+    sourceId: String = "MarchingSquares-source-id",
+    layerId: String = "MarchingSquares-layer-id",
+    palette: (Double?) -> String
+) {
+    mapboxMap.getStyle { style ->
+        val contours = MarchingSquares.extractContours(
+            grid = grid,
+            isoValue = grid.cells.filterNotNull().average()
+        )
+        val features = mutableListOf<Feature>()
+        contours.forEachIndexed { index, polygon: List<Vertex> ->
+            val points = polygon.map {
+                Point.fromLngLat(it.x, it.y)
+            }
+            val lineString = LineString.fromLngLats(points)
+            val feature = Feature.fromGeometry(lineString)
+            val elevation = (index / contours.size) * 255.0
+            feature.addStringProperty("color", palette(elevation))
+            feature.addNumberProperty("elevation", elevation)
+            feature.addStringProperty("type", "contour-line")
+            features.add(feature)
+        }
+        setupContourLayers(
+            style = style,
+            sourceId = sourceId,
+            layerId = layerId,
+            features = features
+        )
+    }
+}

app/src/main/java/com/icegps/geotools/MarchingSquaresContourGenerator.kt

@@ -0,0 +1,194 @@
+package com.icegps.geotools
+
+import com.mapbox.geojson.Point
+import kotlin.math.floor
+
+/**
+ * @author tabidachinokaze
+ * @date 2025/11/21
+ */
+/**
+ * Marching Squares 等高线生成算法
+ */
+class MarchingSquaresContourGenerator {
+
+    fun generateContours(
+        grid: GridModel,
+        contoursCount: Int = 6,
+    ): List<ContourLine> {
+        val elevations = grid.cells.filterNotNull()
+        if (elevations.isEmpty()) return emptyList()
+        val minElev = elevations.minOrNull() ?: 0.0
+        val maxElev = elevations.maxOrNull() ?: 0.0
+
+        val contourInterval = (maxElev - minElev) / contoursCount
+
+        return generateContours(grid, contourInterval, minElev, maxElev)
+    }
+
+
+    fun generateContours(
+        grid: GridModel,
+        contourInterval: Double,
+        minElevation: Double,
+        maxElevation: Double
+    ): List<ContourLine> {
+        val elevations = grid.cells.filterNotNull()
+        if (elevations.isEmpty()) return emptyList()
+
+        val contours = mutableListOf<ContourLine>()
+
+        // 生成等高线层级
+        var currentElev = (floor(minElevation / contourInterval) + 1) * contourInterval
+        while (currentElev <= maxElevation) {
+            val contourLines = generateContourForLevel(grid, currentElev)
+            if (contourLines.isNotEmpty()) {
+                contours.add(ContourLine(currentElev, contourLines))
+            }
+            currentElev += contourInterval
+        }
+
+        return contours
+    }
+
+    private fun generateContourForLevel(grid: GridModel, level: Double): List<List<Point>> {
+        val segments = mutableListOf<List<Point>>()
+
+        // 遍历所有2x2网格
+        for (y in 0 until grid.rows - 1) {
+            for (x in 0 until grid.cols - 1) {
+                // 获取四个角点
+                val corners = floatArrayOf(
+                    grid.getValue(y, x)?.toFloat() ?: continue,
+                    grid.getValue(y, x + 1)?.toFloat() ?: continue,
+                    grid.getValue(y + 1, x + 1)?.toFloat() ?: continue,
+                    grid.getValue(y + 1, x)?.toFloat() ?: continue
+                )
+
+                // 计算配置索引
+                val configIndex = calculateConfiguration(corners, level.toFloat())
+
+                // 根据配置生成线段
+                val cellSegments = generateSegmentsForConfig(grid, x, y, corners, level.toFloat(), configIndex)
+                segments.addAll(cellSegments)
+            }
+        }
+
+        // 连接线段
+        return connectSegments(segments)
+    }
+
+    private fun calculateConfiguration(corners: FloatArray, level: Float): Int {
+        var config = 0
+        for (i in corners.indices) {
+            if (corners[i] >= level) {
+                config = config or (1 shl i)
+            }
+        }
+        return config
+    }
+
+    private fun generateSegmentsForConfig(
+        grid: GridModel,
+        x: Int,
+        y: Int,
+        corners: FloatArray,
+        level: Float,
+        config: Int
+    ): List<List<Point>> {
+        val segments = mutableListOf<List<Point>>()
+
+        // Marching Squares 查找表
+        when (config) {
+            0, 15 -> { /* 无等高线 */
+            }
+
+            else -> {
+                // 计算边上的交点
+                val intersections = mutableListOf<Point>()
+
+                // 检查每条边
+                for (edge in edges) {
+                    val (edgeStart, edgeEnd) = edge
+                    if (hasIntersection(config, edgeStart, edgeEnd)) {
+                        val point = calculateIntersection(grid, x, y, corners, level, edgeStart, edgeEnd)
+                        intersections.add(point)
+                    }
+                }
+
+                // 根据交点生成线段
+                when (intersections.size) {
+                    2 -> segments.add(listOf(intersections[0], intersections[1]))
+                    4 -> {
+                        // 鞍点情况
+                        segments.add(listOf(intersections[0], intersections[1]))
+                        segments.add(listOf(intersections[2], intersections[3]))
+                    }
+                }
+            }
+        }
+
+        return segments
+    }
+
+    private fun hasIntersection(config: Int, edgeStart: Int, edgeEnd: Int): Boolean {
+        val startBit = (config shr edgeStart) and 1
+        val endBit = (config shr edgeEnd) and 1
+        return startBit != endBit
+    }
+
+    private fun calculateIntersection(
+        grid: GridModel,
+        x: Int,
+        y: Int,
+        corners: FloatArray,
+        level: Float,
+        edgeStart: Int,
+        edgeEnd: Int
+    ): Point {
+        val startValue = corners[edgeStart]
+        val endValue = corners[edgeEnd]
+        val t = (level - startValue) / (endValue - startValue)
+
+        // 边上的坐标
+        val (x1, y1) = getEdgePoint(x, y, edgeStart)
+        val (x2, y2) = getEdgePoint(x, y, edgeEnd)
+
+        val interpX = x1 + (x2 - x1) * t
+        val interpY = y1 + (y2 - y1) * t
+
+        return gridToWorld(grid, interpX, interpY)
+    }
+
+    private fun getEdgePoint(x: Int, y: Int, corner: Int): Pair<Double, Double> {
+        return when (corner) {
+            0 -> Pair(x.toDouble(), y.toDouble())
+            1 -> Pair((x + 1).toDouble(), y.toDouble())
+            2 -> Pair((x + 1).toDouble(), (y + 1).toDouble())
+            3 -> Pair(x.toDouble(), (y + 1).toDouble())
+            else -> Pair(x.toDouble(), y.toDouble())
+        }
+    }
+
+    private fun gridToWorld(grid: GridModel, x: Double, y: Double): Point {
+        val lon = grid.minLon + x * (grid.maxLon - grid.minLon) / (grid.cols - 1)
+        val lat = grid.maxLat - y * (grid.maxLat - grid.minLat) / (grid.rows - 1) // Y轴取反
+        return Point.fromLngLat(lon, lat)
+    }
+
+    private fun connectSegments(segments: List<List<Point>>): List<List<Point>> {
+        // 使用之前的连接算法
+        val connector = ContourConnector()
+        return connector.connectSegments(segments)
+    }
+
+    companion object {
+        // 边定义：0-上, 1-右, 2-下, 3-左
+        private val edges = listOf(
+            Pair(0, 1), // 上边
+            Pair(1, 2), // 右边
+            Pair(2, 3), // 下边
+            Pair(3, 0)  // 左边
+        )
+    }
+}

app/src/main/java/com/icegps/geotools/OverallSlopeTrend.kt

@@ -0,0 +1,529 @@
+package com.icegps.geotools
+
+import com.icegps.math.geometry.Angle
+import com.mapbox.geojson.Feature
+import com.mapbox.geojson.FeatureCollection
+import com.mapbox.geojson.LineString
+import com.mapbox.geojson.Point
+import com.mapbox.geojson.Polygon
+import com.mapbox.maps.MapView
+import com.mapbox.maps.Style
+import com.mapbox.maps.extension.style.expressions.generated.Expression
+import com.mapbox.maps.extension.style.layers.addLayer
+import com.mapbox.maps.extension.style.layers.generated.FillLayer
+import com.mapbox.maps.extension.style.layers.generated.LineLayer
+import com.mapbox.maps.extension.style.layers.properties.generated.LineCap
+import com.mapbox.maps.extension.style.layers.properties.generated.LineJoin
+import com.mapbox.maps.extension.style.sources.addSource
+import com.mapbox.maps.extension.style.sources.generated.geoJsonSource
+import kotlin.math.atan2
+import kotlin.math.cos
+import kotlin.math.sin
+import kotlin.math.sqrt
+
+/**
+ * 计算整体地势趋势（平均坡度方向）
+ * 返回角度（0-360度，0表示北，90表示东）和平均坡度
+ */
+fun GridModel.calculateOverallSlopeTrend(): SlopeTrend {
+    var sumDzDx = 0.0
+    var sumDzDy = 0.0
+    var validCount = 0
+
+    for (r in 1 until rows - 1) {  // 跳过边界
+        for (c in 1 until cols - 1) {
+            val slopeDirection = calculateSlopeDirection(r, c)
+            if (slopeDirection != null) {
+                // 将方向转换为向量分量
+                val rad = Math.toRadians(slopeDirection)
+                sumDzDx += cos(rad)
+                sumDzDy += sin(rad)
+                validCount++
+            }
+        }
+    }
+
+    if (validCount == 0) {
+        return SlopeTrend(0.0, 0.0, 0)
+    }
+
+    // 计算平均向量方向
+    val avgDzDx = sumDzDx / validCount
+    val avgDzDy = sumDzDy / validCount
+
+    var overallDirection = Math.toDegrees(atan2(avgDzDy, avgDzDx))
+    if (overallDirection < 0) overallDirection += 360.0
+
+    // 计算平均坡度强度（向量长度）
+    val slopeStrength = sqrt(avgDzDx * avgDzDx + avgDzDy * avgDzDy)
+
+    return SlopeTrend(overallDirection, slopeStrength, validCount)
+}
+
+/**
+ * 坡度趋势结果
+ */
+data class SlopeTrend(
+    val direction: Double,    // 整体趋势方向（度）
+    val strength: Double,     // 趋势强度（0-1）
+    val sampleCount: Int      // 有效样本数
+) {
+    fun getDirectionName(): String {
+        return when {
+            direction >= 337.5 || direction < 22.5 -> "北"
+            direction >= 22.5 && direction < 67.5 -> "东北"
+            direction >= 67.5 && direction < 112.5 -> "东"
+            direction >= 112.5 && direction < 157.5 -> "东南"
+            direction >= 157.5 && direction < 202.5 -> "南"
+            direction >= 202.5 && direction < 247.5 -> "西南"
+            direction >= 247.5 && direction < 292.5 -> "西"
+            else -> "西北"
+        }
+    }
+
+    override fun toString(): String {
+        return "整体趋势: ${getDirectionName()} (${"%.1f".format(direction)}°), 强度: ${"%.3f".format(strength)}"
+    }
+}
+
+// 在Mapbox中显示整体趋势箭头
+/*fun MapView.displayOverallTrendArrow(
+    grid: GridModel,
+    trendSourceId: String,
+    trendLayerId: String
+) {
+    mapboxMap.getStyle { style ->
+        val trend = grid.calculateOverallSlopeTrend()
+
+        if (trend.sampleCount > 0 && trend.strength > 0.1) { // 只有明显趋势才显示
+            val centerLon = (grid.minLon + grid.maxLon) / 2
+            val centerLat = (grid.minLat + grid.maxLat) / 2
+
+            // 计算箭头长度（基于区域大小和趋势强度）
+            val regionWidth = grid.maxLon - grid.minLon
+            val arrowLength = regionWidth * 0.3 * trend.strength
+
+            val arrowDirectionRad = Math.toRadians(trend.direction)
+            val endLon = centerLon + cos(arrowDirectionRad) * arrowLength
+            val endLat = centerLat + sin(arrowDirectionRad) * arrowLength
+
+            // 创建趋势箭头
+            val arrowLine = LineString.fromLngLats(
+                listOf(
+                    Point.fromLngLat(centerLon, centerLat),
+                    Point.fromLngLat(endLon, endLat)
+                )
+            )
+
+            val arrowFeature = Feature.fromGeometry(arrowLine)
+            arrowFeature.addStringProperty("color", "#0000FF") // 蓝色表示整体趋势
+            arrowFeature.addStringProperty("type", "overall-trend")
+            arrowFeature.addNumberProperty("strength", trend.strength)
+
+            val features = listOf(arrowFeature)
+
+            // 添加到地图（处理图层逻辑...）
+            setupTrendLayer(style, trendSourceId, trendLayerId, features)
+
+            println("显示整体趋势箭头: $trend")
+        }
+    }
+}*/
+
+
+/**
+ * 设置趋势箭头图层
+ */
+private fun MapView.setupTrendLayer(
+    style: Style,
+    trendSourceId: String,
+    trendLayerId: String,
+    features: List<Feature>
+) {
+    // 1. 创建数据源
+    val trendSource = geoJsonSource(trendSourceId) {
+        featureCollection(FeatureCollection.fromFeatures(features))
+    }
+
+    // 2. 先移除图层再移除数据源（正确的顺序）
+    try {
+        style.removeStyleLayer(trendLayerId)
+    } catch (_: Exception) {
+    }
+
+    try {
+        style.removeStyleLayer("$trendLayerId-head")
+    } catch (_: Exception) {
+    }
+
+    if (style.styleSourceExists(trendSourceId)) {
+        style.removeStyleSource(trendSourceId)
+    }
+
+    // 3. 添加数据源
+    style.addSource(trendSource)
+
+    // 4. 添加箭头线图层
+    val lineLayer = LineLayer(trendLayerId, trendSourceId).apply {
+        lineColor(Expression.toColor(Expression.get("color")))
+        lineWidth(4.0)  // 整体趋势箭头更粗
+        lineCap(LineCap.ROUND)
+        lineJoin(LineJoin.ROUND)
+        // 可以根据趋势强度调整透明度
+    }
+    style.addLayer(lineLayer)
+
+    // 5. 添加箭头头部图层（三角形）
+    val headLayer = FillLayer("$trendLayerId-head", trendSourceId).apply {
+        fillColor(Expression.toColor(Expression.get("color")))
+    }
+    style.addLayer(headLayer)
+}
+
+fun MapView.displayControllableArrow(
+    grid: GridModel,
+    trendSourceId: String = "controllable-source-id-0",
+    trendLayerId: String = "controllable-layer-id-0",
+    arrowScale: Double = 0.4,
+    angle: Angle,
+    onHeadArrowChange: (List<Point>) -> Unit
+) {
+    mapboxMap.getStyle { style ->
+        val centerLon = (grid.minLon + grid.maxLon) / 2
+        val centerLat = (grid.minLat + grid.maxLat) / 2
+
+        // 计算箭头长度（基于区域大小和趋势强度）
+        val regionWidth = grid.maxLon - grid.minLon
+        val arrowLength = regionWidth * arrowScale * 1.0
+
+        val arrowDirectionRad = angle.radians
+        val endLon = centerLon + sin(arrowDirectionRad) * arrowLength  // 经度用sin
+        val endLat = centerLat + cos(arrowDirectionRad) * arrowLength  // 纬度用cos
+
+        // 创建趋势箭杆
+        val arrowLine = LineString.fromLngLats(
+            listOf(
+                Point.fromLngLat(centerLon, centerLat),
+                Point.fromLngLat(endLon, endLat)
+            )
+        )
+
+        val arrowFeature = Feature.fromGeometry(arrowLine)
+        arrowFeature.addStringProperty("color", "#0000FF") // 蓝色表示整体趋势
+        arrowFeature.addStringProperty("type", "overall-trend")
+
+        // 创建箭头头部
+        val headSize = arrowLength * 0.2
+        val leftRad = arrowDirectionRad + Math.PI * 0.8
+        val rightRad = arrowDirectionRad - Math.PI * 0.8
+
+        val leftLon = endLon + sin(leftRad) * headSize
+        val leftLat = endLat + cos(leftRad) * headSize
+        val rightLon = endLon + sin(rightRad) * headSize
+        val rightLat = endLat + cos(rightRad) * headSize
+
+        val headRing = listOf(
+            Point.fromLngLat(endLon, endLat),
+            Point.fromLngLat(leftLon, leftLat),
+            Point.fromLngLat(rightLon, rightLat),
+            Point.fromLngLat(endLon, endLat)
+        )
+        onHeadArrowChange(headRing)
+        val headPolygon = Polygon.fromLngLats(listOf(headRing))
+        val headFeature = Feature.fromGeometry(headPolygon)
+        headFeature.addStringProperty("color", "#0000FF")
+        headFeature.addStringProperty("type", "overall-trend")
+
+        val features = listOf(arrowFeature, headFeature)
+
+        // 设置图层
+        setupTrendLayer(style, trendSourceId, trendLayerId, features)
+    }
+}
+
+/**
+ * 完整的显示整体趋势箭头方法
+ */
+fun MapView.displayOverallTrendArrow(
+    grid: GridModel,
+    trendSourceId: String,
+    trendLayerId: String,
+    arrowScale: Double = 0.4, // 箭头大小相对于区域尺寸的比例
+    onArrowCenterChange: (Point) -> Unit,
+    onHeadArrowChange: (List<Point>) -> Unit
+) {
+    mapboxMap.getStyle { style ->
+        val trend = grid.calculateOverallSlopeTrend()
+
+        if (trend.sampleCount > 0 && trend.strength > 0.1) { // 只有明显趋势才显示
+            val centerLon = (grid.minLon + grid.maxLon) / 2
+            val centerLat = (grid.minLat + grid.maxLat) / 2
+            onArrowCenterChange(Point.fromLngLat(centerLon, centerLat))
+
+            // 计算箭头长度（基于区域大小和趋势强度）
+            val regionWidth = grid.maxLon - grid.minLon
+            val arrowLength = regionWidth * arrowScale * trend.strength
+
+            val arrowDirectionRad = Math.toRadians(trend.direction)
+            val endLon = centerLon + cos(arrowDirectionRad) * arrowLength
+            val endLat = centerLat + sin(arrowDirectionRad) * arrowLength
+
+            // 创建趋势箭杆
+            val arrowLine = LineString.fromLngLats(
+                listOf(
+                    Point.fromLngLat(centerLon, centerLat),
+                    Point.fromLngLat(endLon, endLat)
+                )
+            )
+
+            val arrowFeature = Feature.fromGeometry(arrowLine)
+            arrowFeature.addStringProperty("color", "#0000FF") // 蓝色表示整体趋势
+            arrowFeature.addStringProperty("type", "overall-trend")
+            arrowFeature.addNumberProperty("strength", trend.strength)
+
+            // 创建箭头头部
+            val headSize = arrowLength * 0.2
+            val leftRad = arrowDirectionRad + Math.PI * 0.8
+            val rightRad = arrowDirectionRad - Math.PI * 0.8
+
+            val leftLon = endLon + cos(leftRad) * headSize
+            val leftLat = endLat + sin(leftRad) * headSize
+            val rightLon = endLon + cos(rightRad) * headSize
+            val rightLat = endLat + sin(rightRad) * headSize
+
+            val headRing = listOf(
+                Point.fromLngLat(endLon, endLat),
+                Point.fromLngLat(leftLon, leftLat),
+                Point.fromLngLat(rightLon, rightLat),
+                Point.fromLngLat(endLon, endLat)
+            )
+            onHeadArrowChange(headRing)
+            val headPolygon = Polygon.fromLngLats(listOf(headRing))
+            val headFeature = Feature.fromGeometry(headPolygon)
+            headFeature.addStringProperty("color", "#0000FF")
+            headFeature.addStringProperty("type", "overall-trend")
+            headFeature.addNumberProperty("strength", trend.strength)
+
+            val features = listOf(arrowFeature, headFeature)
+
+            // 设置图层
+            setupTrendLayer(style, trendSourceId, trendLayerId, features)
+
+            println("显示整体趋势箭头: $trend")
+        } else {
+            // 移除趋势箭头（如果趋势不明显）
+            try {
+                style.removeStyleLayer(trendLayerId)
+            } catch (_: Exception) {
+            }
+
+            try {
+                style.removeStyleLayer("$trendLayerId-head")
+            } catch (_: Exception) {
+            }
+
+            if (style.styleSourceExists(trendSourceId)) {
+                style.removeStyleSource(trendSourceId)
+            }
+
+            println("趋势不明显，已移除趋势箭头")
+        }
+    }
+}
+
+/**
+ * 增强版：显示多个趋势分析结果
+ */
+fun MapView.displayMultipleTrendAnalyses(
+    grid: GridModel,
+    trendSourceId: String = "multiple-trends",
+    trendLayerId: String = "multiple-trends-layer"
+) {
+    mapboxMap.getStyle { style ->
+        val features = mutableListOf<Feature>()
+        val centerLon = (grid.minLon + grid.maxLon) / 2
+        val centerLat = (grid.minLat + grid.maxLat) / 2
+        val regionWidth = grid.maxLon - grid.minLon
+
+        // 1. 平均坡度向量趋势（蓝色）
+        val trend1 = grid.calculateOverallSlopeTrend()
+        if (trend1.strength > 0.1) {
+            addTrendArrow(
+                features, centerLon, centerLat, regionWidth, trend1,
+                "#0000FF", "平均向量", 0.0
+            )
+        }
+
+        // 2. 主导方向趋势（绿色）- 稍微偏移位置
+        val trend2 = grid.calculateDominantSlopeTrend()
+        if (trend2.strength > 0.1) {
+            addTrendArrow(
+                features, centerLon + regionWidth * 0.1, centerLat,
+                regionWidth, trend2, "#00FF00", "主导方向", 1.0
+            )
+        }
+
+        // 3. 梯度分析趋势（红色）- 稍微偏移位置
+        val trend3 = grid.calculateGradientBasedTrend()
+        if (trend3.strength > 0.1) {
+            addTrendArrow(
+                features, centerLon - regionWidth * 0.1, centerLat,
+                regionWidth, trend3, "#FF0000", "梯度分析", 2.0
+            )
+        }
+
+        if (features.isNotEmpty()) {
+            setupTrendLayer(style, trendSourceId, trendLayerId, features)
+            println("显示多个趋势分析结果")
+        } else {
+            // 清理图层
+            try {
+                style.removeStyleLayer(trendLayerId)
+            } catch (_: Exception) {
+            }
+            try {
+                style.removeStyleLayer("$trendLayerId-head")
+            } catch (_: Exception) {
+            }
+            if (style.styleSourceExists(trendSourceId)) {
+                style.removeStyleSource(trendSourceId)
+            }
+        }
+    }
+}
+
+/**
+ * 添加单个趋势箭头到要素列表
+ */
+private fun addTrendArrow(
+    features: MutableList<Feature>,
+    centerLon: Double, centerLat: Double,
+    regionWidth: Double,
+    trend: SlopeTrend,
+    color: String,
+    type: String,
+    offset: Double
+) {
+    val arrowLength = regionWidth * 0.3 * trend.strength
+    val arrowDirectionRad = Math.toRadians(trend.direction)
+
+    val endLon = centerLon + cos(arrowDirectionRad) * arrowLength
+    val endLat = centerLat + sin(arrowDirectionRad) * arrowLength
+
+    // 箭杆
+    val arrowLine = LineString.fromLngLats(
+        listOf(
+            Point.fromLngLat(centerLon, centerLat),
+            Point.fromLngLat(endLon, endLat)
+        )
+    )
+
+    val arrowFeature = Feature.fromGeometry(arrowLine)
+    arrowFeature.addStringProperty("color", color)
+    arrowFeature.addStringProperty("type", type)
+    arrowFeature.addNumberProperty("strength", trend.strength)
+    arrowFeature.addNumberProperty("offset", offset)
+    features.add(arrowFeature)
+
+    // 箭头头部
+    val headSize = arrowLength * 0.15
+    val leftRad = arrowDirectionRad + Math.PI * 0.8
+    val rightRad = arrowDirectionRad - Math.PI * 0.8
+
+    val leftLon = endLon + cos(leftRad) * headSize
+    val leftLat = endLat + sin(leftRad) * headSize
+    val rightLon = endLon + cos(rightRad) * headSize
+    val rightLat = endLat + sin(rightRad) * headSize
+
+    val headRing = listOf(
+        Point.fromLngLat(endLon, endLat),
+        Point.fromLngLat(leftLon, leftLat),
+        Point.fromLngLat(rightLon, rightLat),
+        Point.fromLngLat(endLon, endLat)
+    )
+    val headPolygon = Polygon.fromLngLats(listOf(headRing))
+    val headFeature = Feature.fromGeometry(headPolygon)
+    headFeature.addStringProperty("color", color)
+    headFeature.addStringProperty("type", type)
+    headFeature.addNumberProperty("strength", trend.strength)
+    headFeature.addNumberProperty("offset", offset)
+    features.add(headFeature)
+}
+
+/**
+ * 基于高程梯度分析整体趋势
+ */
+fun GridModel.calculateGradientBasedTrend(): SlopeTrend {
+    var sumGradientX = 0.0
+    var sumGradientY = 0.0
+    var validCount = 0
+
+    for (r in 1 until rows - 1) {
+        for (c in 1 until cols - 1) {
+            val center = getValue(r, c) ?: continue
+            val east = getValue(r, c + 1)
+            val west = getValue(r, c - 1)
+            val north = getValue(r - 1, c)
+            val south = getValue(r + 1, c)
+
+            if (east != null && west != null && north != null && south != null) {
+                // 计算梯度
+                val dzdx = (east - west) / (2 * cellSizeMeters)
+                val dzdy = (north - south) / (2 * cellSizeMeters)  // 注意Y轴方向
+
+                sumGradientX += dzdx
+                sumGradientY += dzdy
+                validCount++
+            }
+        }
+    }
+
+    if (validCount == 0) {
+        return SlopeTrend(0.0, 0.0, 0)
+    }
+
+    val avgGradientX = sumGradientX / validCount
+    val avgGradientY = sumGradientY / validCount
+
+    // 梯度方向是上坡方向，转换为下坡方向
+    var gradientDirection = Math.toDegrees(atan2(avgGradientY, avgGradientX))
+    if (gradientDirection < 0) gradientDirection += 360.0
+
+    val downhillDirection = (gradientDirection + 180.0) % 360.0
+    val strength = sqrt(avgGradientX * avgGradientX + avgGradientY * avgGradientY)
+
+    return SlopeTrend(downhillDirection, strength, validCount)
+}
+
+/**
+ * 使用方向统计计算主导趋势
+ */
+fun GridModel.calculateDominantSlopeTrend(): SlopeTrend {
+    // 将360度分为8个方向区间
+    val directionBins = IntArray(8)  // 0:北, 1:东北, 2:东, 3:东南, 4:南, 5:西南, 6:西, 7:西北
+    var totalCount = 0
+
+    for (r in 1 until rows - 1) {
+        for (c in 1 until cols - 1) {
+            val direction = calculateSlopeDirection(r, c)
+            if (direction != null) {
+                val binIndex = ((direction + 22.5) % 360 / 45).toInt()
+                directionBins[binIndex]++
+                totalCount++
+            }
+        }
+    }
+
+    if (totalCount == 0) {
+        return SlopeTrend(0.0, 0.0, 0)
+    }
+
+    // 找到最多的方向区间
+    val maxBinIndex = directionBins.indices.maxByOrNull { directionBins[it] } ?: 0
+    val dominantDirection = maxBinIndex * 45.0  // 转换为角度
+
+    // 计算主导方向的强度（该方向的比例）
+    val strength = directionBins[maxBinIndex].toDouble() / totalCount
+
+    return SlopeTrend(dominantDirection, strength, totalCount)
+}

app/src/main/java/com/icegps/geotools/RasterUtils.kt

@@ -40,7 +40,7 @@ object RasterUtils {
      * @return Pair(FloatArray(values row-major), Bitmap visualization)
      */
     fun <T : IPoint> rasterizeDelaunay(
-        delaunator: Delaunator<T>,
+        delaunator: IDelaunator<T>,
         minX: Double,
         minY: Double,
         maxX: Double,

app/src/main/java/com/icegps/geotools/RayCastingAlgorithm.kt

@@ -0,0 +1,57 @@
+package com.icegps.geotools
+
+import com.mapbox.geojson.Point
+
+/**
+ * @author tabidachinokaze
+ * @date 2025/11/20
+ */
+/**
+ * 射线法判断点是否在多边形内
+ */
+object RayCastingAlgorithm {
+
+    /**
+     * 使用射线法判断点是否在多边形内
+     * @param point 测试点
+     * @param polygon 多边形顶点列表
+     * @return true如果在多边形内
+     */
+    fun isPointInPolygon(point: Point, polygon: List<Point>): Boolean {
+        if (polygon.size < 3) return false
+
+        val x = point.longitude()
+        val y = point.latitude()
+        var inside = false
+
+        var j = polygon.size - 1
+        for (i in polygon.indices) {
+            val xi = polygon[i].longitude()
+            val yi = polygon[i].latitude()
+            val xj = polygon[j].longitude()
+            val yj = polygon[j].latitude()
+
+            val intersect = ((yi > y) != (yj > y)) &&
+                    (x < (xj - xi) * (y - yi) / (yj - yi) + xi)
+
+            if (intersect) inside = !inside
+            j = i
+        }
+
+        return inside
+    }
+
+    /**
+     * 判断点是否在环内（闭合多边形）
+     */
+    fun isPointInRing(point: Point, ring: List<Point>): Boolean {
+        // 确保环是闭合的（首尾点相同）
+        val closedRing = if (ring.first() != ring.last()) {
+            ring + ring.first()
+        } else {
+            ring
+        }
+
+        return isPointInPolygon(point, closedRing)
+    }
+}

app/src/main/java/com/icegps/geotools/SharedHttpClient.kt

@@ -0,0 +1,41 @@
+package com.icegps.geotools
+
+import io.ktor.client.HttpClient
+import io.ktor.client.engine.cio.CIO
+import io.ktor.client.plugins.HttpTimeout
+import io.ktor.client.plugins.contentnegotiation.ContentNegotiation
+import io.ktor.client.plugins.logging.LogLevel
+import io.ktor.client.plugins.logging.Logger
+import io.ktor.client.plugins.logging.Logging
+import io.ktor.client.plugins.logging.SIMPLE
+import io.ktor.http.ContentType
+import io.ktor.serialization.kotlinx.json.json
+import kotlinx.serialization.json.Json
+
+/**
+ * @author tabidachinokaze
+ * @date 2025/11/20
+ */
+fun SharedHttpClient(json: Json): HttpClient {
+    return HttpClient(CIO) {
+        install(ContentNegotiation) {
+            json(
+                json = json,
+                contentType = ContentType.Text.Html
+            )
+            json(
+                json = json,
+                contentType = ContentType.Application.Json
+            )
+        }
+        install(Logging) {
+            this.level = LogLevel.ALL
+            this.logger = Logger.SIMPLE
+        }
+        install(HttpTimeout) {
+            requestTimeoutMillis = 1000 * 60 * 10
+            connectTimeoutMillis = 1000 * 60 * 5
+            socketTimeoutMillis = 1000 * 60 * 10
+        }
+    }
+}

app/src/main/java/com/icegps/geotools/SharedJson.kt

@@ -0,0 +1,16 @@
+package com.icegps.geotools
+
+import kotlinx.serialization.json.Json
+import kotlinx.serialization.modules.SerializersModule
+
+/**
+ * @author tabidachinokaze
+ * @date 2025/11/20
+ */
+fun SharedJson(): Json {
+    return Json {
+        ignoreUnknownKeys = true
+        serializersModule = SerializersModule {
+        }
+    }
+}

app/src/main/java/com/icegps/geotools/SimpleContourGenerator.kt

@@ -0,0 +1,83 @@
+package com.icegps.geotools
+
+import com.mapbox.geojson.Point
+import kotlin.math.ceil
+
+/**
+ * @author tabidachinokaze
+ * @date 2025/11/21
+ */
+/**
+ * 简化的等高线生成（基于网格边界检测）
+ */
+class SimpleContourGenerator {
+
+    /**
+     * 生成简化的等高线
+     */
+    fun generateSimpleContours(
+        grid: GridModel,
+        contourInterval: Double
+    ): List<ContourLine> {
+        val contours = mutableListOf<ContourLine>()
+        val minElev = grid.cells.filterNotNull().minOrNull() ?: 0.0
+        val maxElev = grid.cells.filterNotNull().maxOrNull() ?: 100.0
+
+        var currentElev = ceil(minElev / contourInterval) * contourInterval
+        while (currentElev <= maxElev) {
+            val lines = generateContourLines(grid, currentElev)
+            if (lines.isNotEmpty()) {
+                contours.add(ContourLine(currentElev, lines))
+            }
+            currentElev += contourInterval
+        }
+
+        return contours
+    }
+
+    private fun generateContourLines(grid: GridModel, elevation: Double): List<List<Point>> {
+        val lines = mutableListOf<List<Point>>()
+
+        // 检查水平边界
+        for (r in 0 until grid.rows) {
+            for (c in 0 until grid.cols - 1) {
+                val left = grid.getValue(r, c) ?: continue
+                val right = grid.getValue(r, c + 1) ?: continue
+
+                if ((left - elevation) * (right - elevation) < 0) {
+                    val t = (elevation - left) / (right - left)
+                    val x = c + t
+                    val y = r.toDouble()
+                    val point = gridToWorld(grid, x, y)
+                    // 简化：每个交点作为单独的点
+                    lines.add(listOf(point))
+                }
+            }
+        }
+
+        // 检查垂直边界
+        for (r in 0 until grid.rows - 1) {
+            for (c in 0 until grid.cols) {
+                val top = grid.getValue(r, c) ?: continue
+                val bottom = grid.getValue(r + 1, c) ?: continue
+
+                if ((top - elevation) * (bottom - elevation) < 0) {
+                    val t = (elevation - top) / (bottom - top)
+                    val x = c.toDouble()
+                    val y = r + t
+                    val point = gridToWorld(grid, x, y)
+                    lines.add(listOf(point))
+                }
+            }
+        }
+
+        return lines
+    }
+
+    private fun gridToWorld(grid: GridModel, gridX: Double, gridY: Double): Point {
+        val lon = grid.minLon + gridX * (grid.maxLon - grid.minLon) / grid.cols
+        // 修正：Y轴方向取反，因为row=0对应北边（maxLat），row增大向南（minLat）
+        val lat = grid.maxLat - gridY * (grid.maxLat - grid.minLat) / grid.rows
+        return Point.fromLngLat(lon, lat)
+    }
+}

app/src/main/java/com/icegps/geotools/TerrainAnalyzer.kt

@@ -0,0 +1,20 @@
+package com.icegps.geotools
+
+/**
+ * 栅格数据模型
+ */
+data class GridDEM(
+    val cols: Int,
+    val rows: Int,
+    val xllCorner: Int,
+    val yllCorner: Int,
+    val cellSize: Double
+)
+
+/**
+ * @author tabidachinokaze
+ * @date 2025/11/19
+ */
+class TerrainAnalyzer {
+}
+

app/src/main/java/com/icegps/geotools/VolcanoGenerator.kt

@@ -0,0 +1,294 @@
+package com.icegps.geotools
+
+import com.icegps.math.geometry.Vector3D
+import kotlin.math.PI
+import kotlin.math.abs
+import kotlin.math.cos
+import kotlin.math.exp
+import kotlin.math.floor
+import kotlin.math.max
+import kotlin.math.min
+import kotlin.math.sin
+import kotlin.math.sqrt
+
+/**
+ * @author tabidachinokaze
+ * @date 2025/11/19
+ */
+class VolcanoGenerator {
+
+    // 基础火山锥形
+    fun generateVolcanoCone(
+        radius: Int,
+        height: Double = 50.0,
+        resolution: Int = 100
+    ): List<Vector3D> {
+        val points = mutableListOf<Vector3D>()
+
+        for (i in 0 until resolution) {
+            val angle = 2.0 * PI * i / resolution
+            for (r in 0..radius step 2) {
+                val currentHeight = height * (1.0 - r.toDouble() / radius)
+                // 添加一些噪声使表面更自然
+                val noise = perlinNoise(r * cos(angle), r * sin(angle), 0.2) * 2.0
+
+                val x = r * cos(angle)
+                val z = r * sin(angle)
+                val y = max(0.0, currentHeight + noise)
+
+                points.add(Vector3D(x, y, z))
+            }
+        }
+
+        return points
+    }
+
+    // 复杂火山地形（包含火山口）
+    fun generateComplexVolcano(
+        baseRadius: Int = 50,
+        craterRadius: Int = 15,
+        maxHeight: Double = 80.0,
+        resolution: Int = 120
+    ): List<Vector3D> {
+        val points = mutableListOf<Vector3D>()
+
+        for (i in 0 until resolution) {
+            val angle = 2.0 * PI * i / resolution
+
+            for (r in 0..baseRadius step 2) {
+                val normalizedR = r.toDouble() / baseRadius
+
+                // 火山剖面函数
+                val height = when {
+                    r <= craterRadius -> {
+                        // 火山口内部
+                        val craterDepth = maxHeight * 0.3
+                        craterDepth * (1.0 - normalizedR * 2)
+                    }
+
+                    else -> {
+                        // 火山锥外部
+                        val coneHeight = maxHeight * exp(-normalizedR * 2)
+                        coneHeight
+                    }
+                }
+
+                // 添加侵蚀效果和噪声
+                val erosion = perlinNoise(r * cos(angle) * 0.1, r * sin(angle) * 0.1, 0.1) * 5.0
+                val detailNoise = perlinNoise(r * cos(angle), r * sin(angle), 0.05) * 3.0
+
+                val x = r * cos(angle)
+                val z = r * sin(angle)
+                val y = max(0.0, height + erosion + detailNoise)
+
+                points.add(Vector3D(x, y, z))
+            }
+        }
+
+        return points
+    }
+
+    // 火山群地形
+    fun generateVolcanoField(
+        width: Int,
+        depth: Int,
+        scale: Double = 2.0,
+        volcanoCount: Int = 3
+    ): List<Vector3D> {
+        val points = mutableListOf<Vector3D>()
+        val volcanoes = generateVolcanoPositions(volcanoCount, width, depth, scale)
+
+        for (x in 0 until width) {
+            for (z in 0 until depth) {
+                val xPos = x * scale
+                val zPos = z * scale
+                var height = 0.0
+
+                // 叠加多个火山
+                for (volcano in volcanoes) {
+                    val dx = xPos - volcano.x
+                    val dz = zPos - volcano.z
+                    val distance = sqrt(dx * dx + dz * dz)
+
+                    if (distance <= volcano.radius) {
+                        val volcanoHeight = calculateVolcanoHeight(distance, volcano)
+                        height += volcanoHeight
+                    }
+                }
+
+                // 添加基础地形噪声
+                val baseNoise = perlinNoise(xPos * 0.02, zPos * 0.02, 0.1) * 3.0
+                points.add(Vector3D(xPos, height + baseNoise, zPos))
+            }
+        }
+
+        return points
+    }
+
+    // 活跃火山（带有熔岩流）
+    fun generateActiveVolcano(
+        baseRadius: Int = 60,
+        maxHeight: Double = 100.0,
+        lavaFlowCount: Int = 4
+    ): List<Vector3D> {
+        val points = mutableListOf<Vector3D>()
+        val lavaFlows = generateLavaFlowPaths(lavaFlowCount, baseRadius)
+
+        for (i in 0 until 360 step 2) {
+            val angle = Math.toRadians(i.toDouble())
+
+            for (r in 0..baseRadius step 2) {
+                val normalizedR = r.toDouble() / baseRadius
+
+                // 基础火山形状
+                var height = when {
+                    r <= 20 -> {
+                        // 深火山口
+                        maxHeight * 0.7 * (1.0 - normalizedR)
+                    }
+
+                    else -> {
+                        // 陡峭的火山锥
+                        maxHeight * exp(-normalizedR * 1.5)
+                    }
+                }
+
+                // 添加熔岩流效果
+                val lavaEffect = calculateLavaFlowEffect(r, angle, lavaFlows, maxHeight)
+                height += lavaEffect
+
+                // 地形细节
+                val erosion = perlinNoise(r * cos(angle) * 0.15, r * sin(angle) * 0.15, 0.08) * 8.0
+                val detail = perlinNoise(r * cos(angle), r * sin(angle), 0.03) * 4.0
+
+                val x = r * cos(angle)
+                val z = r * sin(angle)
+                val y = max(0.0, height + erosion + detail)
+
+                points.add(Vector3D(x, y, z))
+            }
+        }
+
+        return points
+    }
+
+    // 海底火山
+    fun generateUnderwaterVolcano(
+        baseRadius: Int = 40,
+        maxHeight: Double = 60.0,
+        waterDepth: Double = 30.0
+    ): List<Vector3D> {
+        val points = mutableListOf<Vector3D>()
+
+        for (i in 0 until 360 step 3) {
+            val angle = Math.toRadians(i.toDouble())
+
+            for (r in 0..baseRadius step 2) {
+                val normalizedR = r.toDouble() / baseRadius
+
+                // 水下火山形状（更平缓）
+                var height = maxHeight * exp(-normalizedR * 3.0) - waterDepth
+
+                // 添加水下侵蚀效果
+                val underwaterErosion = perlinNoise(r * cos(angle) * 0.2, r * sin(angle) * 0.2, 0.1) * 6.0
+                val sediment = perlinNoise(r * cos(angle) * 0.05, r * sin(angle) * 0.05, 0.02) * 2.0
+
+                val x = r * cos(angle)
+                val z = r * sin(angle)
+                val y = height + underwaterErosion + sediment
+
+                points.add(Vector3D(x, y, z))
+            }
+        }
+
+        return points
+    }
+
+    // 辅助函数
+    private data class VolcanoInfo(val x: Double, val z: Double, val radius: Double, val maxHeight: Double)
+
+    private data class LavaFlow(val startAngle: Double, val width: Double, val intensity: Double)
+
+    private fun generateVolcanoPositions(count: Int, width: Int, depth: Int, scale: Double): List<VolcanoInfo> {
+        return List(count) {
+            VolcanoInfo(
+                x = (width * 0.2 + width * 0.6 * Math.random()) * scale,
+                z = (depth * 0.2 + depth * 0.6 * Math.random()) * scale,
+                radius = (20.0 + Math.random() * 30.0),
+                maxHeight = (40.0 + Math.random() * 60.0)
+            )
+        }
+    }
+
+    private fun generateLavaFlowPaths(count: Int, baseRadius: Int): List<LavaFlow> {
+        return List(count) {
+            LavaFlow(
+                startAngle = Math.random() * 2 * PI,
+                width = 0.2 + Math.random() * 0.3,
+                intensity = 5.0 + Math.random() * 15.0
+            )
+        }
+    }
+
+    private fun calculateVolcanoHeight(distance: Double, volcano: VolcanoInfo): Double {
+        val normalizedDistance = distance / volcano.radius
+        return when {
+            normalizedDistance < 0.3 -> {
+                // 火山口区域
+                volcano.maxHeight * 0.8 * (1.0 - normalizedDistance * 2)
+            }
+
+            else -> {
+                // 火山锥区域
+                volcano.maxHeight * exp(-normalizedDistance * 2)
+            }
+        }
+    }
+
+    private fun calculateLavaFlowEffect(r: Int, angle: Double, lavaFlows: List<LavaFlow>, maxHeight: Double): Double {
+        var effect = 0.0
+        val normalizedR = r.toDouble() / 60.0
+
+        for (flow in lavaFlows) {
+            var angleDiff = abs(angle - flow.startAngle)
+            angleDiff = min(angleDiff, 2 * PI - angleDiff)
+
+            if (angleDiff < flow.width) {
+                // 熔岩流效果随距离衰减
+                val flowEffect = flow.intensity * exp(-normalizedR * 3) * (1.0 - angleDiff / flow.width)
+                effect += flowEffect
+            }
+        }
+
+        return effect
+    }
+
+    private fun perlinNoise(x: Double, z: Double, frequency: Double = 0.1): Double {
+        val x0 = floor(x * frequency)
+        val z0 = floor(z * frequency)
+        val x1 = x0 + 1
+        val z1 = z0 + 1
+
+        fun grad(ix: Int, iz: Int): Double {
+            val random = sin(ix * 12.9898 + iz * 78.233) * 43758.5453
+            return (random % 1.0) * 2 - 1
+        }
+
+        fun interpolate(a: Double, b: Double, w: Double): Double {
+            return a + (b - a) * (w * w * (3 - 2 * w))
+        }
+
+        val g00 = grad(x0.toInt(), z0.toInt())
+        val g10 = grad(x1.toInt(), z0.toInt())
+        val g01 = grad(x0.toInt(), z1.toInt())
+        val g11 = grad(x1.toInt(), z1.toInt())
+
+        val tx = x * frequency - x0
+        val tz = z * frequency - z0
+
+        val n0 = interpolate(g00, g10, tx)
+        val n1 = interpolate(g01, g11, tx)
+
+        return interpolate(n0, n1, tz)
+    }
+}

app/src/main/java/com/icegps/geotools/api/LookupResponse.kt

@@ -0,0 +1,11 @@
+package com.icegps.geotools.api
+
+import com.icegps.geotools.model.GeoPoint
+import kotlinx.serialization.SerialName
+import kotlinx.serialization.Serializable
+
+@Serializable
+data class LookupResponse(
+    @SerialName("results")
+    val results: List<GeoPoint>
+)

app/src/main/java/com/icegps/geotools/api/OpenElevationApi.kt

@@ -0,0 +1,36 @@
+package com.icegps.geotools.api
+
+import com.icegps.geotools.model.IPoint
+import com.icegps.geotools.model.latitude
+import com.icegps.geotools.model.longitude
+import io.ktor.client.HttpClient
+import io.ktor.client.call.body
+import io.ktor.client.request.get
+import io.ktor.client.request.parameter
+import io.ktor.http.appendPathSegments
+import io.ktor.http.path
+
+/**
+ * @author tabidachinokaze
+ * @date 2025/11/20
+ */
+interface OpenElevationApi {
+    suspend fun lookup(values: List<IPoint>): LookupResponse
+}
+
+class OpenElevation(
+    private val client: HttpClient
+) : OpenElevationApi {
+    private val baseUrl: String = "https://api.open-elevation.com/api/v1/"
+
+    // curl 'https://api.open-elevation.com/api/v1/lookup?locations=10,10|20,20|41.161758,-8.583933'
+    override suspend fun lookup(values: List<IPoint>): LookupResponse {
+        val response = client.get(baseUrl) {
+            url {
+                appendPathSegments("lookup")
+                parameter("locations", values.joinToString("|") { "${it.latitude},${it.longitude}" })
+            }
+        }
+        return response.body()
+    }
+}

app/src/main/java/com/icegps/geotools/ktx/Context.kt

@@ -0,0 +1,12 @@
+package com.icegps.geotools.ktx
+
+import android.content.Context
+import android.widget.Toast
+
+/**
+ * @author tabidachinokaze
+ * @date 2025/11/20
+ */
+fun Context.toast(text: String, duration: Int = Toast.LENGTH_SHORT) {
+    Toast.makeText(this, text, duration).show()
+}

app/src/main/java/com/icegps/geotools/model/DPoint.kt

@@ -1,5 +1,9 @@
 package com.icegps.geotools.model
 
+interface IGeoPoint : IPoint {
+    var z: Double
+}
+
 /**
  * @author tabidachinokaze
  * @date 2025/11/5
@@ -7,5 +11,5 @@ package com.icegps.geotools.model
 data class DPoint(
     override var x: Double,
     override var y: Double,
-    var z: Double
-) : IPoint
+    override var z: Double
+) : IGeoPoint

app/src/main/java/com/icegps/geotools/model/GeoPoint.kt

@@ -0,0 +1,34 @@
+package com.icegps.geotools.model
+
+import kotlinx.serialization.SerialName
+import kotlinx.serialization.Serializable
+
+/**
+ * @author tabidachinokaze
+ * @date 2025/11/20
+ */
+@Serializable
+data class GeoPoint(
+    @SerialName("latitude")
+    var latitude: Double,
+    @SerialName("longitude")
+    var longitude: Double,
+    @SerialName("elevation")
+    var elevation: Double
+) : IGeoPoint {
+    override var z: Double
+        get() = elevation
+        set(value) {
+            elevation = value
+        }
+    override var x: Double
+        get() = longitude
+        set(value) {
+            longitude = value
+        }
+    override var y: Double
+        get() = latitude
+        set(value) {
+            latitude = value
+        }
+}

app/src/main/res/layout/activity_main.xml

@@ -1,10 +1,116 @@
 <?xml version="1.0" encoding="utf-8"?>
-<androidx.constraintlayout.widget.ConstraintLayout xmlns:android="http://schemas.android.com/apk/res/android"
-    xmlns:app="http://schemas.android.com/apk/res-auto"
+<LinearLayout xmlns:android="http://schemas.android.com/apk/res/android"
     xmlns:tools="http://schemas.android.com/tools"
     android:id="@+id/main"
     android:layout_width="match_parent"
     android:layout_height="match_parent"
+    android:orientation="horizontal"
     tools:context=".MainActivity">
 
-</androidx.constraintlayout.widget.ConstraintLayout>
+
+    <LinearLayout
+        android:layout_width="0dp"
+        android:layout_height="match_parent"
+        android:layout_weight="1"
+        android:orientation="vertical">
+
+        <LinearLayout
+            android:layout_width="wrap_content"
+            android:layout_height="wrap_content"
+            android:gravity="center_vertical">
+
+            <TextView
+                android:layout_width="wrap_content"
+                android:layout_height="wrap_content"
+                android:text="栅格大小：" />
+
+            <com.google.android.material.slider.Slider
+                android:id="@+id/slider"
+                android:layout_width="match_parent"
+                android:layout_height="wrap_content"
+                android:value="50"
+                android:valueFrom="50"
+                android:valueTo="500" />
+        </LinearLayout>
+
+        <LinearLayout
+            android:layout_width="wrap_content"
+            android:layout_height="wrap_content"
+            android:gravity="center_vertical">
+
+            <TextView
+                android:layout_width="wrap_content"
+                android:layout_height="wrap_content"
+                android:text="角度(0~360)：" />
+
+            <com.google.android.material.slider.Slider
+                android:id="@+id/angle_slider"
+                android:layout_width="match_parent"
+                android:layout_height="wrap_content"
+                android:value="0"
+                android:valueFrom="0"
+                android:valueTo="360" />
+        </LinearLayout>
+
+        <LinearLayout
+            android:layout_width="wrap_content"
+            android:layout_height="wrap_content"
+            android:gravity="center_vertical">
+
+            <TextView
+                android:layout_width="wrap_content"
+                android:layout_height="wrap_content"
+                android:text="坡度(%)：" />
+
+            <com.google.android.material.slider.Slider
+                android:id="@+id/slope_percentage_slider"
+                android:layout_width="match_parent"
+                android:layout_height="wrap_content"
+                android:value="50"
+                android:valueFrom="0"
+                android:valueTo="100" />
+        </LinearLayout>
+
+        <LinearLayout
+            android:layout_width="wrap_content"
+            android:layout_height="wrap_content"
+            android:gravity="center_vertical">
+
+            <TextView
+                android:layout_width="wrap_content"
+                android:layout_height="wrap_content"
+                android:text="基准高度(m)：" />
+
+            <com.google.android.material.slider.Slider
+                android:id="@+id/base_height_offset_slider"
+                android:layout_width="match_parent"
+                android:layout_height="wrap_content"
+                android:value="0"
+                android:valueFrom="-100"
+                android:valueTo="100" />
+        </LinearLayout>
+
+        <LinearLayout
+            android:layout_width="match_parent"
+            android:layout_height="wrap_content">
+
+            <Button
+                android:id="@+id/button_clear"
+                android:layout_width="wrap_content"
+                android:layout_height="wrap_content"
+                android:text="清除点" />
+        </LinearLayout>
+
+        <TextView
+            android:id="@+id/design_cut"
+            android:layout_width="match_parent"
+            android:layout_height="wrap_content"
+            android:text="设计土方量：" />
+    </LinearLayout>
+
+    <com.mapbox.maps.MapView
+        android:id="@+id/map_view"
+        android:layout_width="0dp"
+        android:layout_height="match_parent"
+        android:layout_weight="3" />
+</LinearLayout>

app/src/test/java/com/icegps/geotools/FindContoursTest.kt

@@ -0,0 +1,57 @@
+package com.icegps.geotools
+
+import com.icegps.common.helper.GeoHelper
+import com.icegps.geotools.model.DPoint
+import com.icegps.geotools.model.IGeoPoint
+import kotlinx.coroutines.runBlocking
+import kotlinx.coroutines.withTimeout
+import org.junit.Test
+import org.openrndr.math.Vector2
+import org.openrndr.shape.LineSegment
+import org.openrndr.shape.Rectangle
+import org.openrndr.shape.ShapeContour
+import kotlin.math.max
+import kotlin.math.min
+
+/**
+ * @author tabidachinokaze
+ * @date 2025/11/21
+ */
+class FindContoursTest {
+    @Test
+    fun testFindContours() = runBlocking {
+        withTimeout(10000) {
+            val geoHelper = GeoHelper.getSharedInstance()
+            initGeoHelper()
+            val points = coordinateGenerate1()
+            val delaunator: IDelaunator<IGeoPoint> = Delaunator(
+                points = points.map {
+                    val wgs84 = geoHelper.enuToWGS84Object(GeoHelper.ENU(it.x, it.y, it.z))
+                    DPoint(wgs84.lon, wgs84.lat, wgs84.hgt)
+                }
+            )
+            val cellSizeMeters = 2.0
+            val gridModel = triangulationToGrid(
+                delaunator = delaunator,
+                cellSizeMeters = cellSizeMeters,
+            )
+
+            val rectangle = Rectangle(
+                x = gridModel.minLon,
+                y = gridModel.minLat,
+                width = gridModel.maxLon - gridModel.minLon,
+                height = gridModel.maxLat - gridModel.minLat
+            )
+
+            val contours = findContours(
+                grid = gridModel,
+            )
+            // org.openrndr.extra.marchingsquares.findContours()
+            println("contours size = ${contours.size}")
+            contours.forEach {
+                println(it)
+            }
+        }
+        Unit
+    }
+}

app/src/test/java/com/icegps/geotools/MathTest.kt

@@ -0,0 +1,41 @@
+package com.icegps.geotools
+
+import org.junit.Test
+import kotlin.math.log2
+import kotlin.math.pow
+
+/**
+ * @author tabidachinokaze
+ * @date 2025/11/21
+ */
+class MathTest {
+    @Test
+    fun testPow() {
+        (0..100).map { index ->
+            index.toDouble().pow(2) * 0.0008
+        }.forEach {
+            println(it)
+        }
+    }
+
+    @Test
+    fun testLog() {
+        (0..100).map { index ->
+            println(log2(index.toDouble()))
+        }
+    }
+
+    @Test
+    fun testLog2() {
+        (0..100).map { x ->
+            (0..100).map { y ->
+                println(log2(x * y * 1.0) * 2.0)
+            }
+        }
+    }
+
+    @Test
+    fun testEarthwork() {
+
+    }
+}

app/src/test/java/com/icegps/geotools/OpenElevationTest.kt

@@ -0,0 +1,24 @@
+package com.icegps.geotools
+
+import com.icegps.geotools.api.OpenElevation
+import com.icegps.geotools.api.OpenElevationApi
+import com.icegps.geotools.model.GeoPoint
+import kotlinx.coroutines.runBlocking
+import org.junit.Test
+
+/**
+ * @author tabidachinokaze
+ * @date 2025/11/20
+ */
+class OpenElevationTest {
+    private val openElevation: OpenElevationApi = OpenElevation(
+        client = SharedHttpClient(SharedJson())
+    )
+
+    @Test
+    fun testOpenElevation() = runBlocking {
+        val response = openElevation.lookup(listOf(GeoPoint(41.161758, -8.583933, 0.0)))
+        println(response)
+        Unit
+    }
+}

build.gradle

@@ -4,4 +4,5 @@ plugins {
     alias(libs.plugins.kotlin.android) apply false
     alias(libs.plugins.kotlin.jvm) apply false
     alias(libs.plugins.android.library) apply false
+    alias(libs.plugins.kotlin.serialization) apply false
 }

delaunator/src/main/java/com/icegps/geotools/model/IPoint.kt

@@ -4,3 +4,6 @@ interface IPoint {
     var x: Double
     var y: Double
 }
+
+val IPoint.longitude: Double get() = x
+val IPoint.latitude: Double get() = y

gradle/libs.versions.toml

@@ -9,7 +9,8 @@ appcompat = "1.7.1"
 material = "1.12.0"
 activity = "1.11.0"
 constraintlayout = "2.2.1"
-
+kotlinx-serialization = "1.7.0"
+ktor = "2.3.6"
 [libraries]
 androidx-core-ktx = { group = "androidx.core", name = "core-ktx", version.ref = "coreKtx" }
 junit = { group = "junit", name = "junit", version.ref = "junit" }
@@ -20,10 +21,16 @@ material = { group = "com.google.android.material", name = "material", version.r
 androidx-activity = { group = "androidx.activity", name = "activity", version.ref = "activity" }
 androidx-constraintlayout = { group = "androidx.constraintlayout", name = "constraintlayout", version.ref = "constraintlayout" }
 kotlin-test = { module = "org.jetbrains.kotlin:kotlin-test", version.ref = "kotlin" }
-
+kotlinx-serialization-json = { group = "org.jetbrains.kotlinx", name = "kotlinx-serialization-json", version.ref = "kotlinx-serialization" }
+#ktor
+ktor-client-core = { group = "io.ktor", name = "ktor-client-core", version.ref = "ktor" }
+ktor-client-cio = { group = "io.ktor", name = "ktor-client-cio", version.ref = "ktor" }
+ktor-client-logging = { group = "io.ktor", name = "ktor-client-logging", version.ref = "ktor" }
+ktor-client-content-negotiation = { group = "io.ktor", name = "ktor-client-content-negotiation", version.ref = "ktor" }
+ktor-serialization-kotlinx-json = { group = "io.ktor", name = "ktor-serialization-kotlinx-json", version.ref = "ktor" }
 [plugins]
 android-application = { id = "com.android.application", version.ref = "agp" }
 kotlin-android = { id = "org.jetbrains.kotlin.android", version.ref = "kotlin" }
 kotlin-jvm = { id = "org.jetbrains.kotlin.jvm", version.ref = "kotlin" }
 android-library = { id = "com.android.library", version.ref = "agp" }
-
+kotlin-serialization = { id = "org.jetbrains.kotlin.plugin.serialization", version.ref = "kotlin" }