From a096f6ee0869a55f6a0cbdc08f85e9fdc03055c2 Mon Sep 17 00:00:00 2001 From: tabidachinokaze Date: Fri, 21 Nov 2025 20:12:41 +0800 Subject: [PATCH] =?UTF-8?q?feat:=20=E5=AE=9E=E7=8E=B0=E6=89=8B=E5=8A=BF?= =?UTF-8?q?=E8=B0=83=E6=95=B4=E5=9D=A1=E5=BA=A6?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- app/build.gradle | 19 +- .../com/icegps/geotools/ContourConnector.kt | 402 ++++++++ .../com/icegps/geotools/ContourGenerator.kt | 327 ++++++ .../icegps/geotools/DraggableTrendArrow.kt | 314 ++++++ .../com/icegps/geotools/EarthworkManager.kt | 937 ++++++++++++++++++ .../java/com/icegps/geotools/FindContours.kt | 319 ++++++ .../java/com/icegps/geotools/GeoJsonUtils.kt | 4 +- .../main/java/com/icegps/geotools/GridCell.kt | 496 ++++++++- .../geotools/HeightmapVolcanoGenerator.kt | 375 +++++++ .../java/com/icegps/geotools/MainActivity.kt | 470 ++++++++- .../com/icegps/geotools/MarchingSquares.kt | 321 ++++++ .../MarchingSquaresContourGenerator.kt | 194 ++++ .../com/icegps/geotools/OverallSlopeTrend.kt | 529 ++++++++++ .../java/com/icegps/geotools/RasterUtils.kt | 2 +- .../icegps/geotools/RayCastingAlgorithm.kt | 57 ++ .../com/icegps/geotools/SharedHttpClient.kt | 41 + .../java/com/icegps/geotools/SharedJson.kt | 16 + .../icegps/geotools/SimpleContourGenerator.kt | 83 ++ .../com/icegps/geotools/TerrainAnalyzer.kt | 20 + .../com/icegps/geotools/VolcanoGenerator.kt | 294 ++++++ .../com/icegps/geotools/api/LookupResponse.kt | 11 + .../icegps/geotools/api/OpenElevationApi.kt | 36 + .../java/com/icegps/geotools/ktx/Context.kt | 12 + .../java/com/icegps/geotools/model/DPoint.kt | 8 +- .../com/icegps/geotools/model/GeoPoint.kt | 34 + app/src/main/res/layout/activity_main.xml | 112 ++- .../com/icegps/geotools/FindContoursTest.kt | 57 ++ .../test/java/com/icegps/geotools/MathTest.kt | 41 + .../com/icegps/geotools/OpenElevationTest.kt | 24 + build.gradle | 1 + .../java/com/icegps/geotools/model/IPoint.kt | 5 +- gradle/libs.versions.toml | 13 +- 32 files changed, 5476 insertions(+), 98 deletions(-) create mode 100644 app/src/main/java/com/icegps/geotools/ContourConnector.kt create mode 100644 app/src/main/java/com/icegps/geotools/ContourGenerator.kt create mode 100644 app/src/main/java/com/icegps/geotools/DraggableTrendArrow.kt create mode 100644 app/src/main/java/com/icegps/geotools/EarthworkManager.kt create mode 100644 app/src/main/java/com/icegps/geotools/FindContours.kt create mode 100644 app/src/main/java/com/icegps/geotools/HeightmapVolcanoGenerator.kt create mode 100644 app/src/main/java/com/icegps/geotools/MarchingSquares.kt create mode 100644 app/src/main/java/com/icegps/geotools/MarchingSquaresContourGenerator.kt create mode 100644 app/src/main/java/com/icegps/geotools/OverallSlopeTrend.kt create mode 100644 app/src/main/java/com/icegps/geotools/RayCastingAlgorithm.kt create mode 100644 app/src/main/java/com/icegps/geotools/SharedHttpClient.kt create mode 100644 app/src/main/java/com/icegps/geotools/SharedJson.kt create mode 100644 app/src/main/java/com/icegps/geotools/SimpleContourGenerator.kt create mode 100644 app/src/main/java/com/icegps/geotools/TerrainAnalyzer.kt create mode 100644 app/src/main/java/com/icegps/geotools/VolcanoGenerator.kt create mode 100644 app/src/main/java/com/icegps/geotools/api/LookupResponse.kt create mode 100644 app/src/main/java/com/icegps/geotools/api/OpenElevationApi.kt create mode 100644 app/src/main/java/com/icegps/geotools/ktx/Context.kt create mode 100644 app/src/main/java/com/icegps/geotools/model/GeoPoint.kt create mode 100644 app/src/test/java/com/icegps/geotools/FindContoursTest.kt create mode 100644 app/src/test/java/com/icegps/geotools/MathTest.kt create mode 100644 app/src/test/java/com/icegps/geotools/OpenElevationTest.kt diff --git a/app/build.gradle b/app/build.gradle index 16ada42..922ffb5 100644 --- a/app/build.gradle +++ b/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 diff --git a/app/src/main/java/com/icegps/geotools/ContourConnector.kt b/app/src/main/java/com/icegps/geotools/ContourConnector.kt new file mode 100644 index 0000000..87d6f2a --- /dev/null +++ b/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> { + if (segments.isEmpty()) return emptyList() + + val continuousContours = mutableListOf>() + 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>, + usedSegments: BooleanArray, + startIndex: Int + ): List { + val contour = mutableListOf() + 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): 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 { + 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() + + 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, tolerance: Double): List { + if (contour.size <= 2) return contour + + return douglasPeucker(contour, 0, contour.size - 1, tolerance) + } + + /** + * 道格拉斯-普克算法实现 + */ + private fun douglasPeucker( + points: List, + start: Int, + end: Int, + tolerance: Double + ): List { + 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> { + val allLines = mutableListOf>() + + 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> { + 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>() + val intersections = mutableListOf() + + // 检查四条边 + 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() + + 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 条连续线段") + } +} \ No newline at end of file diff --git a/app/src/main/java/com/icegps/geotools/ContourGenerator.kt b/app/src/main/java/com/icegps/geotools/ContourGenerator.kt new file mode 100644 index 0000000..c43c1b9 --- /dev/null +++ b/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 { + val elevations = grid.cells.filterNotNull() + val minElev = minElevation ?: elevations.minOrNull() ?: 0.0 + val maxElev = maxElevation ?: elevations.maxOrNull() ?: 100.0 + + val contours = mutableListOf() + + // 生成等高线层级 + 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> { + val contours = mutableListOf>() + 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 + ): List { + val contour = mutableListOf() + 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 + ): List { + val intersections = mutableListOf() + 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>, // 等高线路径（可能有多个闭合环） + val properties: Map = emptyMap() +) { + fun getLineStrings(): List { + 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() + + // 为每条等高线创建要素 + 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 +) { + 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) +} \ No newline at end of file diff --git a/app/src/main/java/com/icegps/geotools/DraggableTrendArrow.kt b/app/src/main/java/com/icegps/geotools/DraggableTrendArrow.kt new file mode 100644 index 0000000..d5de316 --- /dev/null +++ b/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) { + 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 +} \ No newline at end of file diff --git a/app/src/main/java/com/icegps/geotools/EarthworkManager.kt b/app/src/main/java/com/icegps/geotools/EarthworkManager.kt new file mode 100644 index 0000000..0fa24d6 --- /dev/null +++ b/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, + slopePercentage: Double, + slopeDirection: Double + ): GridModel { + val designCells = Array(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 { + val solutions = mutableListOf() + + // 方案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, + 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 { + 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() + + 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() + 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() + mapboxMap.getStyle { style -> + val features = mutableListOf() + 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, + 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, + 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 { + 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, + 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(val first: A, val second: B, val third: C, val fourth: D) \ No newline at end of file diff --git a/app/src/main/java/com/icegps/geotools/FindContours.kt b/app/src/main/java/com/icegps/geotools/FindContours.kt new file mode 100644 index 0000000..5a0ba9b --- /dev/null +++ b/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 { + val segments = mutableListOf() + val segmentsMap = mutableMapOf>() + + 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() + val contours = mutableListOf() + for (segment in segments) { + if (segment in processedSegments) { + continue + } else { + val collected = mutableListOf() + 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() + 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() + + // 添加调试信息 + 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 +) { + 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) +} \ No newline at end of file diff --git a/app/src/main/java/com/icegps/geotools/GeoJsonUtils.kt b/app/src/main/java/com/icegps/geotools/GeoJsonUtils.kt index 06025c3..fef5cf9 100644 --- a/app/src/main/java/com/icegps/geotools/GeoJsonUtils.kt +++ b/app/src/main/java/com/icegps/geotools/GeoJsonUtils.kt @@ -13,7 +13,7 @@ import com.mapbox.geojson.Polygon */ object GeoJsonUtils { // 生成三角形（Polygon）FeatureCollection - fun trianglesToPolygons(delaunator: Delaunator): FeatureCollection { + fun trianglesToPolygons(delaunator: IDelaunator): FeatureCollection { val features = mutableListOf() val tris = delaunator.triangles // triangles 是按 3 个索引为一组三角形存储 @@ -48,7 +48,7 @@ object GeoJsonUtils { } // 生成边（LineString）FeatureCollection（每条边只输出一次） - fun trianglesToUniqueEdges(delaunator: Delaunator): FeatureCollection { + fun trianglesToUniqueEdges(delaunator: IDelaunator): FeatureCollection { val features = mutableListOf() val seen = HashSet>() val tris = delaunator.triangles diff --git a/app/src/main/java/com/icegps/geotools/GridCell.kt b/app/src/main/java/com/icegps/geotools/GridCell.kt index 2e8e3d7..58a40ff 100644 --- a/app/src/main/java/com/icegps/geotools/GridCell.kt +++ b/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 // 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, + delaunator: IDelaunator, cellSizeMeters: Double = 50.0, // 每个格子的边长（米） maxSidePixels: Int = 5000 // 限制 max rows/cols 防止 OOM（可选） ): GridModel { @@ -158,7 +239,7 @@ fun triangulationToGrid( val cells = Array(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 +) { + fun setRange(range: ClosedFloatingPointRange) { + this.range = range + } + + private val colors: Map + + 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 { + val colorMap = mutableMapOf() + + // 定义关键颜色点 + 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() mapboxMap.getStyle { style -> + val polygonFeatures = mutableListOf() + val arrowFeatures = mutableListOf() + val textFeatures = mutableListOf() - val features = mutableListOf() - - // 计算每格经纬度跨度 - 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) } } diff --git a/app/src/main/java/com/icegps/geotools/HeightmapVolcanoGenerator.kt b/app/src/main/java/com/icegps/geotools/HeightmapVolcanoGenerator.kt new file mode 100644 index 0000000..8521cf5 --- /dev/null +++ b/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 { + val points = mutableListOf() + + 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 { + val points = mutableListOf() + 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 { + val points = mutableListOf() + 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 { + val points = mutableListOf() + 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 { + val points = mutableListOf() + + // 在一条线上生成多个火山 + 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): 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 { + 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 { + 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() +} \ No newline at end of file diff --git a/app/src/main/java/com/icegps/geotools/MainActivity.kt b/app/src/main/java/com/icegps/geotools/MainActivity.kt index 7928a91..701b122 100644 --- a/app/src/main/java/com/icegps/geotools/MainActivity.kt +++ b/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 + private val palette = SimplePalette(0.0..255.0) + private val delaunatorFlow = MutableStateFlow?>(null) + private val clickedPoints = MutableStateFlow>(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()) 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 { + 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 { + /** + * 绕 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 { + val lists = (0..10).map { x -> + (0..10).map { y -> + Vector3D(x, y, x * y) + } + }.flatten() + return lists +} + +fun coordinateGenerate4(): List { + 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 { + 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) { +fun MapView.loadGeoJson(delaunator: IDelaunator) { + 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) { 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) { +fun MapView.loadRaterFromTin(delaunator: IDelaunator) { // add ImageSource delaunator.points val bbox = RasterUtils.boundingBox(delaunator.points) @@ -290,10 +698,10 @@ fun MapView.loadRaterFromTin(delaunator: Delaunator) { valueGetter = { it.z } ) val bboxCoordinates: List> = 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, + delaunator: IDelaunator, testSourceId: String, testLayerId: String ) { @@ -460,7 +868,7 @@ fun MapView.showVoronoiAsPolygons( * - pixelsPerDegree: 控制生成的栅格分辨率（每经度多少像素），或直接给 width/height */ fun MapView.rasterizeDelaunayToMap( - delaunator: Delaunator, + delaunator: IDelaunator, testSourceId: String, testLayerId: String, pixelsPerDegree: Double = 400.0 // 可调：值越大图片越精细、越大 diff --git a/app/src/main/java/com/icegps/geotools/MarchingSquares.kt b/app/src/main/java/com/icegps/geotools/MarchingSquares.kt new file mode 100644 index 0000000..4f65a10 --- /dev/null +++ b/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 + +/** + * Marching Squares 算法实现 + * 用于从栅格数据中提取等值线 + */ +object MarchingSquares { + /** + * 从栅格数据中提取等值线 + * + * @param grid 栅格数据模型 + * @param isoValue 等值线对应的数值 + * @return 等值线列表，每个等值线是一个多边形顶点列表 + */ + fun extractContours(grid: GridModel, isoValue: Double): List { + val contours = mutableListOf() + + // 遍历所有网格单元 (注意边界: 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): List { + if (segments.isEmpty()) return emptyList() + + val connectedContours = mutableListOf() + val used = BooleanArray(segments.size) { false } + + for (i in segments.indices) { + if (!used[i] && segments[i].size >= 2) { + val currentContour = mutableListOf() + 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() + contours.forEachIndexed { index, polygon: List -> + 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 + ) + } +} \ No newline at end of file diff --git a/app/src/main/java/com/icegps/geotools/MarchingSquaresContourGenerator.kt b/app/src/main/java/com/icegps/geotools/MarchingSquaresContourGenerator.kt new file mode 100644 index 0000000..35ed19a --- /dev/null +++ b/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 { + 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 { + val elevations = grid.cells.filterNotNull() + if (elevations.isEmpty()) return emptyList() + + val contours = mutableListOf() + + // 生成等高线层级 + 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> { + val segments = mutableListOf>() + + // 遍历所有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> { + val segments = mutableListOf>() + + // Marching Squares 查找表 + when (config) { + 0, 15 -> { /* 无等高线 */ + } + + else -> { + // 计算边上的交点 + val intersections = mutableListOf() + + // 检查每条边 + 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 { + 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> { + // 使用之前的连接算法 + 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) // 左边 + ) + } +} diff --git a/app/src/main/java/com/icegps/geotools/OverallSlopeTrend.kt b/app/src/main/java/com/icegps/geotools/OverallSlopeTrend.kt new file mode 100644 index 0000000..8bfccb8 --- /dev/null +++ b/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 +) { + // 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) -> 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) -> 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() + 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, + 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) +} diff --git a/app/src/main/java/com/icegps/geotools/RasterUtils.kt b/app/src/main/java/com/icegps/geotools/RasterUtils.kt index 5ff8e77..f0d18fb 100644 --- a/app/src/main/java/com/icegps/geotools/RasterUtils.kt +++ b/app/src/main/java/com/icegps/geotools/RasterUtils.kt @@ -40,7 +40,7 @@ object RasterUtils { * @return Pair(FloatArray(values row-major), Bitmap visualization) */ fun rasterizeDelaunay( - delaunator: Delaunator, + delaunator: IDelaunator, minX: Double, minY: Double, maxX: Double, diff --git a/app/src/main/java/com/icegps/geotools/RayCastingAlgorithm.kt b/app/src/main/java/com/icegps/geotools/RayCastingAlgorithm.kt new file mode 100644 index 0000000..a287b10 --- /dev/null +++ b/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): 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): Boolean { + // 确保环是闭合的（首尾点相同） + val closedRing = if (ring.first() != ring.last()) { + ring + ring.first() + } else { + ring + } + + return isPointInPolygon(point, closedRing) + } +} \ No newline at end of file diff --git a/app/src/main/java/com/icegps/geotools/SharedHttpClient.kt b/app/src/main/java/com/icegps/geotools/SharedHttpClient.kt new file mode 100644 index 0000000..191ffbd --- /dev/null +++ b/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 + } + } +} \ No newline at end of file diff --git a/app/src/main/java/com/icegps/geotools/SharedJson.kt b/app/src/main/java/com/icegps/geotools/SharedJson.kt new file mode 100644 index 0000000..8a23d2d --- /dev/null +++ b/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 { + } + } +} \ No newline at end of file diff --git a/app/src/main/java/com/icegps/geotools/SimpleContourGenerator.kt b/app/src/main/java/com/icegps/geotools/SimpleContourGenerator.kt new file mode 100644 index 0000000..332e63b --- /dev/null +++ b/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 { + val contours = mutableListOf() + 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> { + val lines = mutableListOf>() + + // 检查水平边界 + 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) + } +} \ No newline at end of file diff --git a/app/src/main/java/com/icegps/geotools/TerrainAnalyzer.kt b/app/src/main/java/com/icegps/geotools/TerrainAnalyzer.kt new file mode 100644 index 0000000..1204367 --- /dev/null +++ b/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 { +} + diff --git a/app/src/main/java/com/icegps/geotools/VolcanoGenerator.kt b/app/src/main/java/com/icegps/geotools/VolcanoGenerator.kt new file mode 100644 index 0000000..888ec0f --- /dev/null +++ b/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 { + val points = mutableListOf() + + 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 { + val points = mutableListOf() + + 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 { + val points = mutableListOf() + 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 { + val points = mutableListOf() + 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 { + val points = mutableListOf() + + 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 { + 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 { + 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, 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) + } +} \ No newline at end of file diff --git a/app/src/main/java/com/icegps/geotools/api/LookupResponse.kt b/app/src/main/java/com/icegps/geotools/api/LookupResponse.kt new file mode 100644 index 0000000..4ccf875 --- /dev/null +++ b/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 +) \ No newline at end of file diff --git a/app/src/main/java/com/icegps/geotools/api/OpenElevationApi.kt b/app/src/main/java/com/icegps/geotools/api/OpenElevationApi.kt new file mode 100644 index 0000000..49aad20 --- /dev/null +++ b/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): 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): LookupResponse { + val response = client.get(baseUrl) { + url { + appendPathSegments("lookup") + parameter("locations", values.joinToString("|") { "${it.latitude},${it.longitude}" }) + } + } + return response.body() + } +} \ No newline at end of file diff --git a/app/src/main/java/com/icegps/geotools/ktx/Context.kt b/app/src/main/java/com/icegps/geotools/ktx/Context.kt new file mode 100644 index 0000000..2493ec4 --- /dev/null +++ b/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() +} \ No newline at end of file diff --git a/app/src/main/java/com/icegps/geotools/model/DPoint.kt b/app/src/main/java/com/icegps/geotools/model/DPoint.kt index 0be50ee..3efcfe6 100644 --- a/app/src/main/java/com/icegps/geotools/model/DPoint.kt +++ b/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 diff --git a/app/src/main/java/com/icegps/geotools/model/GeoPoint.kt b/app/src/main/java/com/icegps/geotools/model/GeoPoint.kt new file mode 100644 index 0000000..3a47e68 --- /dev/null +++ b/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 + } +} \ No newline at end of file diff --git a/app/src/main/res/layout/activity_main.xml b/app/src/main/res/layout/activity_main.xml index 9affce0..99062cf 100644 --- a/app/src/main/res/layout/activity_main.xml +++ b/app/src/main/res/layout/activity_main.xml @@ -1,10 +1,116 @@ - - \ No newline at end of file + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +