package de.arxsilex.maptools.converter.core.elevation

import de.arxsilex.maptools.converter.core.ElevationOutputFormat
import de.arxsilex.maptools.converter.core.EARTH_MEAN_RADIUS_KM
import de.arxsilex.maptools.converter.core.GeoBounds
import de.arxsilex.maptools.converter.core.GeoCircleRegion
import de.arxsilex.maptools.converter.core.GeoRectangleRegion
import de.arxsilex.maptools.converter.core.GeoSelectionRegion
import de.arxsilex.maptools.converter.core.InputFileInspector
import de.arxsilex.maptools.converter.core.InputFormat
import de.arxsilex.maptools.converter.core.unionOrNull
import java.awt.image.RenderedImage
import java.io.BufferedOutputStream
import java.io.Closeable
import java.io.DataOutputStream
import java.io.RandomAccessFile
import java.nio.ByteBuffer
import java.nio.ByteOrder
import java.nio.MappedByteBuffer
import java.nio.channels.FileChannel
import java.nio.charset.StandardCharsets
import java.nio.file.AtomicMoveNotSupportedException
import java.nio.file.Files
import java.nio.file.Path
import java.nio.file.StandardCopyOption
import java.nio.file.StandardOpenOption
import java.util.Comparator
import java.util.Locale
import java.util.concurrent.Callable
import java.util.concurrent.ExecutionException
import java.util.concurrent.ExecutorCompletionService
import java.util.concurrent.Executors
import kotlin.io.path.fileSize
import kotlin.io.path.name
import kotlin.math.ceil
import kotlin.math.acos
import kotlin.math.cos
import kotlin.math.floor
import kotlin.math.max
import kotlin.math.PI
import kotlin.math.roundToInt
import kotlin.math.sin
import org.apache.sis.coverage.grid.PixelInCell
import org.apache.sis.geometry.DirectPosition2D
import org.apache.sis.referencing.CRS
import org.apache.sis.referencing.CommonCRS
import org.apache.sis.storage.Aggregate
import org.apache.sis.storage.DataStore
import org.apache.sis.storage.DataStores
import org.apache.sis.storage.GridCoverageResource
import org.opengis.referencing.operation.MathTransform

data class ElevationConversionRequest(
    val inputs: List<Path>,
    val output: Path,
    val outputFormat: ElevationOutputFormat,
    val areas: List<GeoSelectionRegion>,
    val executionOptions: ElevationExecutionOptions = ElevationExecutionPolicy.detect(),
)

data class ElevationConversionReport(
    val outputs: List<Path>,
    val sampledCells: Long,
    val populatedCells: Long,
)

private data class GeoTiffStripResult(
    val processedCells: Long,
    val populatedCells: Long,
)

/** Converts independent DEM sources into one de-duplicated output grid or a set of standard tiles. */
object ElevationConverter {
    fun convert(
        request: ElevationConversionRequest,
        onProgress: (completed: Long, total: Long) -> Unit = { _, _ -> },
        isCancelled: () -> Boolean = { false },
    ): ElevationConversionReport {
        require(request.inputs.isNotEmpty()) { "Mindestens eine Höhendatenquelle ist erforderlich." }
        require(request.areas.isNotEmpty()) { "Mindestens ein Auswahlbereich ist erforderlich." }
        require(!Files.exists(request.output)) { "Das Ausgabeziel existiert bereits: ${request.output}" }
        val sources = request.inputs.map(::openSource).sortedBy(ElevationSource::cellSizeDegrees)
        try {
            val sourceBounds = sources.map(ElevationSource::bounds).unionOrNull()
                ?: error("Die Höhendaten besitzen keine geografischen Grenzen.")
            val areaBounds = request.areas.map(GeoSelectionRegion::bounds).unionOrNull()
                ?: error("Die Auswahl besitzt keine geografischen Grenzen.")
            val outputBounds = intersect(sourceBounds, areaBounds)
                ?: error("Keine Höhendatenquelle überschneidet den gewählten Exportbereich.")
            val sampler = MergedElevationSampler(sources, request.areas)
            return when (request.outputFormat) {
                ElevationOutputFormat.GEOTIFF -> writeGeoTiff(
                    request.output,
                    outputBounds,
                    sampler,
                    sources.first().cellSizeDegrees,
                    request.executionOptions,
                    onProgress,
                    isCancelled,
                )
                ElevationOutputFormat.ESRI_ASCII_GRID -> writeAsciiGrid(
                    request.output,
                    outputBounds,
                    sampler,
                    sources.first().cellSizeDegrees,
                    onProgress,
                    isCancelled,
                )
                ElevationOutputFormat.SRTM_HGT -> writeHgtTiles(
                    request.output,
                    outputBounds,
                    sampler,
                    sources.first().cellSizeDegrees,
                    onProgress,
                    isCancelled,
                )
                ElevationOutputFormat.DTED -> writeDtedTiles(
                    request.output,
                    outputBounds,
                    sampler,
                    sources.first().cellSizeDegrees,
                    onProgress,
                    isCancelled,
                )
            }
        } finally {
            sources.forEach { runCatching { it.close() } }
        }
    }

    private fun writeAsciiGrid(
        output: Path,
        bounds: GeoBounds,
        sampler: MergedElevationSampler,
        requestedCellSize: Double,
        onProgress: (Long, Long) -> Unit,
        isCancelled: () -> Boolean,
    ): ElevationConversionReport {
        val grid = OutputGrid.create(bounds, requestedCellSize)
        val temporary = output.partial()
        Files.createDirectories(output.toAbsolutePath().parent)
        var populated = 0L
        try {
            Files.newBufferedWriter(temporary, StandardCharsets.US_ASCII).use { writer ->
                writer.appendLine("ncols ${grid.width}")
                writer.appendLine("nrows ${grid.height}")
                writer.appendLine("xllcorner ${formatDecimal(grid.bounds.west)}")
                writer.appendLine("yllcorner ${formatDecimal(grid.bounds.south)}")
                writer.appendLine("cellsize ${formatDecimal(grid.cellSize)}")
                writer.appendLine("NODATA_value ${NO_DATA.toInt()}")
                for (row in 0 until grid.height) {
                    checkNotCancelled(isCancelled)
                    val latitude = grid.latitudeAt(row)
                    for (column in 0 until grid.width) {
                        if (column > 0) writer.append(' ')
                        val value = sampler.sample(grid.longitudeAt(column), latitude)
                        if (value == null) {
                            writer.append(NO_DATA.toInt().toString())
                        } else {
                            populated++
                            writer.append(formatElevation(value))
                        }
                    }
                    writer.newLine()
                    onProgress((row + 1L) * grid.width, grid.cellCount)
                }
            }
            moveAtomically(temporary, output)
            return ElevationConversionReport(listOf(output), grid.cellCount, populated)
        } finally {
            Files.deleteIfExists(temporary)
        }
    }

    private fun writeGeoTiff(
        output: Path,
        bounds: GeoBounds,
        sampler: MergedElevationSampler,
        requestedCellSize: Double,
        executionOptions: ElevationExecutionOptions,
        onProgress: (Long, Long) -> Unit,
        isCancelled: () -> Boolean,
    ): ElevationConversionReport {
        val grid = OutputGrid.create(bounds, requestedCellSize)
        val temporary = output.partial()
        Files.createDirectories(output.toAbsolutePath().parent)
        try {
            val populated = writeStreamingGeoTiff(
                temporary,
                grid,
                sampler,
                executionOptions,
                onProgress,
                isCancelled,
            )
            require(Files.isRegularFile(temporary) && Files.size(temporary) > 0L) {
                "Der GeoTIFF-Writer hat keine Datei erzeugt."
            }
            moveAtomically(temporary, output)
            return ElevationConversionReport(listOf(output), grid.cellCount, populated)
        } finally {
            Files.deleteIfExists(temporary)
        }
    }

    /** Writes an uncompressed, strip-oriented WGS84 Float32 GeoTIFF without retaining the raster in RAM. */
    private fun writeStreamingGeoTiff(
        output: Path,
        grid: OutputGrid,
        sampler: MergedElevationSampler,
        executionOptions: ElevationExecutionOptions,
        onProgress: (Long, Long) -> Unit,
        isCancelled: () -> Boolean,
    ): Long {
        val rowBytes = grid.width.toLong() * Float.SIZE_BYTES
        require(rowBytes in 1..MAX_GEOTIFF_STRIP_BYTES.toLong()) {
            "Das GeoTIFF ist mit ${grid.width} Spalten zu breit. Bitte wähle einen schmaleren Bereich."
        }
        val targetStripBytes = executionOptions.bufferBytesPerWorker.coerceAtMost(MAX_GEOTIFF_STRIP_BYTES)
        val rowsPerStrip = (targetStripBytes / rowBytes.toInt()).coerceIn(1, grid.height)
        val stripCount = ((grid.height + rowsPerStrip - 1L) / rowsPerStrip).toInt()
        val entryCount = 15
        var extraCursor = align(8L + 2L + entryCount * 12L + 4L, 8L)
        fun allocate(bytes: Long, alignment: Long): Long {
            extraCursor = align(extraCursor, alignment)
            val offset = extraCursor
            extraCursor += bytes
            return offset
        }

        val pixelScaleOffset = allocate(3L * Double.SIZE_BYTES, 8L)
        val tiePointOffset = allocate(6L * Double.SIZE_BYTES, 8L)
        val geoKeys = shortArrayOf(
            1, 1, 0, 4,
            1024, 0, 1, 2,
            1025, 0, 1, 1,
            2048, 0, 1, 4326,
            2054, 0, 1, 9102,
        )
        val geoKeyOffset = allocate(geoKeys.size.toLong() * Short.SIZE_BYTES, 2L)
        val noDataBytes = "${NO_DATA.toInt()}\u0000".toByteArray(StandardCharsets.US_ASCII)
        val noDataOffset = allocate(noDataBytes.size.toLong(), 1L)
        val stripOffsetsArrayOffset = if (stripCount > 1) allocate(stripCount * 4L, 4L) else 0L
        val stripCountsArrayOffset = if (stripCount > 1) allocate(stripCount * 4L, 4L) else 0L
        val dataOffset = align(extraCursor, 4L)
        val stripByteCounts = LongArray(stripCount) { index ->
            val firstRow = index * rowsPerStrip
            minOf(rowsPerStrip, grid.height - firstRow).toLong() * rowBytes
        }
        val stripOffsets = LongArray(stripCount)
        var nextOffset = dataOffset
        for (index in 0 until stripCount) {
            stripOffsets[index] = nextOffset
            nextOffset += stripByteCounts[index]
        }
        require(nextOffset <= UINT32_MAX) {
            "Das GeoTIFF würde die 4-GiB-Grenze des klassischen TIFF-Formats überschreiten."
        }

        RandomAccessFile(output.toFile(), "rw").use { file ->
            file.setLength(0L)
            file.write(byteArrayOf('I'.code.toByte(), 'I'.code.toByte()))
            file.writeShortLe(42)
            file.writeIntLe(8L)
            file.seek(8L)
            file.writeShortLe(entryCount)
            file.writeIfdEntry(256, TIFF_LONG, 1, grid.width.toLong())
            file.writeIfdEntry(257, TIFF_LONG, 1, grid.height.toLong())
            file.writeIfdEntry(258, TIFF_SHORT, 1, 32)
            file.writeIfdEntry(259, TIFF_SHORT, 1, 1)
            file.writeIfdEntry(262, TIFF_SHORT, 1, 1)
            file.writeIfdEntry(273, TIFF_LONG, stripCount.toLong(), if (stripCount == 1) stripOffsets[0] else stripOffsetsArrayOffset)
            file.writeIfdEntry(277, TIFF_SHORT, 1, 1)
            file.writeIfdEntry(278, TIFF_LONG, 1, rowsPerStrip.toLong())
            file.writeIfdEntry(279, TIFF_LONG, stripCount.toLong(), if (stripCount == 1) stripByteCounts[0] else stripCountsArrayOffset)
            file.writeIfdEntry(284, TIFF_SHORT, 1, 1)
            file.writeIfdEntry(339, TIFF_SHORT, 1, 3)
            file.writeIfdEntry(33550, TIFF_DOUBLE, 3, pixelScaleOffset)
            file.writeIfdEntry(33922, TIFF_DOUBLE, 6, tiePointOffset)
            file.writeIfdEntry(34735, TIFF_SHORT, geoKeys.size.toLong(), geoKeyOffset)
            file.writeIfdEntry(42113, TIFF_ASCII, noDataBytes.size.toLong(), noDataOffset)
            file.writeIntLe(0)

            file.seek(pixelScaleOffset)
            file.writeDoubleLe(grid.cellSize)
            file.writeDoubleLe(grid.cellSize)
            file.writeDoubleLe(0.0)
            file.seek(tiePointOffset)
            listOf(0.0, 0.0, 0.0, grid.bounds.west, grid.bounds.north, 0.0).forEach { value ->
                file.writeDoubleLe(value)
            }
            file.seek(geoKeyOffset)
            geoKeys.forEach { file.writeShortLe(it.toInt()) }
            file.seek(noDataOffset)
            file.write(noDataBytes)
            if (stripCount > 1) {
                file.seek(stripOffsetsArrayOffset)
                stripOffsets.forEach { value -> file.writeIntLe(value) }
                file.seek(stripCountsArrayOffset)
                stripByteCounts.forEach { value -> file.writeIntLe(value) }
            }
            file.setLength(nextOffset)
        }

        val executor = Executors.newFixedThreadPool(executionOptions.workerCount) { runnable ->
            Thread(runnable, "maptools-elevation-geotiff").apply { isDaemon = true }
        }
        val completion = ExecutorCompletionService<GeoTiffStripResult>(executor)
        try {
            FileChannel.open(output, StandardOpenOption.WRITE).use { channel ->
                for (stripIndex in 0 until stripCount) {
                    completion.submit(
                        Callable {
                            writeGeoTiffStrip(
                                channel = channel,
                                fileOffset = stripOffsets[stripIndex],
                                firstRow = stripIndex * rowsPerStrip,
                                rowCount = (stripByteCounts[stripIndex] / rowBytes).toInt(),
                                grid = grid,
                                sampler = sampler,
                                isCancelled = isCancelled,
                            )
                        },
                    )
                }
                var populated = 0L
                var completedCells = 0L
                repeat(stripCount) {
                    val result = try {
                        completion.take().get()
                    } catch (exception: ExecutionException) {
                        throw exception.cause ?: exception
                    }
                    populated += result.populatedCells
                    completedCells += result.processedCells
                    onProgress(completedCells.coerceAtMost(grid.cellCount), grid.cellCount)
                }
                channel.force(false)
                return populated
            }
        } finally {
            executor.shutdownNow()
        }
    }

    private fun writeGeoTiffStrip(
        channel: FileChannel,
        fileOffset: Long,
        firstRow: Int,
        rowCount: Int,
        grid: OutputGrid,
        sampler: MergedElevationSampler,
        isCancelled: () -> Boolean,
    ): GeoTiffStripResult {
        val buffer = ByteBuffer.allocate(Math.multiplyExact(Math.multiplyExact(rowCount, grid.width), Float.SIZE_BYTES))
            .order(ByteOrder.LITTLE_ENDIAN)
        var populated = 0L
        for (row in firstRow until firstRow + rowCount) {
            checkNotCancelled(isCancelled)
            val latitude = grid.latitudeAt(row)
            val selectedIntervals = sampler.selectedLongitudeIntervals(latitude)
            var intervalIndex = 0
            for (column in 0 until grid.width) {
                val longitude = grid.longitudeAt(column)
                while (
                    intervalIndex < selectedIntervals.size &&
                    longitude > selectedIntervals[intervalIndex].endInclusive
                ) {
                    intervalIndex++
                }
                val selected = intervalIndex < selectedIntervals.size &&
                    longitude >= selectedIntervals[intervalIndex].start
                val value = if (selected) sampler.sampleSelected(longitude, latitude) else null
                if (value != null) populated++
                buffer.putFloat(value?.toFloat() ?: NO_DATA)
            }
        }
        buffer.flip()
        var writeOffset = fileOffset
        while (buffer.hasRemaining()) {
            checkNotCancelled(isCancelled)
            val written = channel.write(buffer, writeOffset)
            check(written > 0) { "Der GeoTIFF-Writer konnte den Ausgabepuffer nicht schreiben." }
            writeOffset += written
        }
        return GeoTiffStripResult(
            processedCells = rowCount.toLong() * grid.width,
            populatedCells = populated,
        )
    }

    private fun writeHgtTiles(
        outputDirectory: Path,
        bounds: GeoBounds,
        sampler: MergedElevationSampler,
        requestedCellSize: Double,
        onProgress: (Long, Long) -> Unit,
        isCancelled: () -> Boolean,
    ): ElevationConversionReport {
        val dimension = if (requestedCellSize <= 1.0 / 2400.0) 3601 else 1201
        return writeDegreeTiles(
            outputDirectory,
            bounds,
            dimension,
            extension = "hgt",
            sampler,
            onProgress,
            isCancelled,
        ) { target, south, west, report ->
            writeHgtTile(target, south, west, dimension, sampler, report, isCancelled)
        }
    }

    private fun writeDtedTiles(
        outputDirectory: Path,
        bounds: GeoBounds,
        sampler: MergedElevationSampler,
        requestedCellSize: Double,
        onProgress: (Long, Long) -> Unit,
        isCancelled: () -> Boolean,
    ): ElevationConversionReport {
        val level = when {
            requestedCellSize <= 1.0 / 2400.0 -> 2
            requestedCellSize <= 1.0 / 240.0 -> 1
            else -> 0
        }
        val dimension = when (level) {
            2 -> 3601
            1 -> 1201
            else -> 121
        }
        return writeDegreeTiles(
            outputDirectory,
            bounds,
            dimension,
            extension = "dt$level",
            sampler,
            onProgress,
            isCancelled,
        ) { target, south, west, report ->
            writeDtedTile(target, south, west, dimension, level, sampler, report, isCancelled)
        }
    }

    private fun writeDegreeTiles(
        outputDirectory: Path,
        bounds: GeoBounds,
        dimension: Int,
        extension: String,
        sampler: MergedElevationSampler,
        onProgress: (Long, Long) -> Unit,
        isCancelled: () -> Boolean,
        writeTile: (Path, Int, Int, (Boolean) -> Unit) -> Unit,
    ): ElevationConversionReport {
        val tiles = buildList {
            for (south in floor(bounds.south).toInt() until ceil(bounds.north).toInt()) {
                for (west in floor(bounds.west).toInt() until ceil(bounds.east).toInt()) {
                    add(south to west)
                }
            }
        }
        require(tiles.isNotEmpty()) { "Der gewählte Bereich enthält keine vollständige Höhendatenkachel." }
        val temporary = outputDirectory.resolveSibling("${outputDirectory.fileName}.partial")
        require(!Files.exists(temporary)) { "Temporäres Ausgabeziel existiert bereits: $temporary" }
        Files.createDirectories(outputDirectory.toAbsolutePath().parent)
        Files.createDirectories(temporary)
        val total = tiles.size.toLong() * dimension * dimension
        var completed = 0L
        var populated = 0L
        val outputs = mutableListOf<Path>()
        try {
            tiles.forEach { (south, west) ->
                checkNotCancelled(isCancelled)
                val name = tileName(south, west, extension)
                val target = temporary.resolve(name)
                writeTile(target, south, west) { hasValue ->
                    completed++
                    if (hasValue) populated++
                    if (completed % dimension == 0L || completed == total) onProgress(completed, total)
                }
                outputs.add(outputDirectory.resolve(name))
            }
            moveAtomically(temporary, outputDirectory)
            return ElevationConversionReport(outputs, total, populated)
        } finally {
            if (Files.exists(temporary)) deleteDirectory(temporary)
        }
    }

    private fun writeHgtTile(
        target: Path,
        south: Int,
        west: Int,
        dimension: Int,
        sampler: MergedElevationSampler,
        report: (Boolean) -> Unit,
        isCancelled: () -> Boolean,
    ) {
        val divisor = (dimension - 1).toDouble()
        DataOutputStream(BufferedOutputStream(Files.newOutputStream(target))).use { output ->
            for (row in 0 until dimension) {
                checkNotCancelled(isCancelled)
                val latitude = south + 1.0 - row / divisor
                for (column in 0 until dimension) {
                    val value = sampler.sample(west + column / divisor, latitude)
                    output.writeShort(toSignedElevation(value).toInt())
                    report(value != null)
                }
            }
        }
    }

    private fun writeDtedTile(
        target: Path,
        south: Int,
        west: Int,
        dimension: Int,
        level: Int,
        sampler: MergedElevationSampler,
        report: (Boolean) -> Unit,
        isCancelled: () -> Boolean,
    ) {
        val intervalTenths = when (level) {
            2 -> 10
            1 -> 30
            else -> 300
        }
        DataOutputStream(BufferedOutputStream(Files.newOutputStream(target))).use { output ->
            output.write(dtedHeader(south, west, dimension, intervalTenths))
            val divisor = (dimension - 1).toDouble()
            for (column in 0 until dimension) {
                checkNotCancelled(isCancelled)
                val record = ByteArray(8 + dimension * 2)
                record[0] = 0xAA.toByte()
                record[1] = (column ushr 16).toByte()
                record[2] = (column ushr 8).toByte()
                record[3] = column.toByte()
                record[4] = (column ushr 8).toByte()
                record[5] = column.toByte()
                val buffer = ByteBuffer.wrap(record).order(ByteOrder.BIG_ENDIAN)
                val longitude = west + column / divisor
                for (row in 0 until dimension) {
                    val latitude = south + row / divisor
                    val value = sampler.sample(longitude, latitude)
                    buffer.putShort(8 + row * 2, toSignedMagnitudeElevation(value))
                    report(value != null)
                }
                output.write(record)
                val checksum = record.sumOf { it.toUByte().toLong() } and 0xFFFF_FFFFL
                output.writeInt(checksum.toInt())
            }
        }
    }

    internal fun openSource(path: Path): ElevationSource {
        val format = InputFileInspector.detectFormat(path.name)
            ?: error("Unbekanntes Höhendatenformat: ${path.fileName}")
        require(format.isElevationData) { "${path.fileName} ist keine Höhendatenquelle." }
        val metadata = ElevationMetadataReader.read(path, format)
        return when (format) {
            InputFormat.ELEVATION_HGT -> HgtElevationSource(path, metadata.bounds, metadata.width)
            InputFormat.ELEVATION_ASC -> AsciiGridElevationSource(path)
            InputFormat.ELEVATION_DTED -> DtedElevationSource(path, metadata.bounds, metadata.width, metadata.height)
            InputFormat.ELEVATION_GEOTIFF -> GeoTiffElevationSource(path, metadata.bounds)
            else -> error("Nicht unterstütztes Höhendatenformat: $format")
        }
    }

    private fun intersect(first: GeoBounds, second: GeoBounds): GeoBounds? {
        val west = maxOf(first.west, second.west)
        val south = maxOf(first.south, second.south)
        val east = minOf(first.east, second.east)
        val north = minOf(first.north, second.north)
        return if (west < east && south < north) GeoBounds(west, south, east, north) else null
    }

    private fun checkNotCancelled(isCancelled: () -> Boolean) {
        if (isCancelled()) throw InterruptedException("Die Höhenkonvertierung wurde abgebrochen.")
    }

    private fun formatDecimal(value: Double): String = String.format(Locale.ROOT, "%.12f", value)
        .trimEnd('0')
        .trimEnd('.')

    private fun formatElevation(value: Double): String {
        val rounded = value.roundToInt()
        return if (kotlin.math.abs(value - rounded) < 0.000_5) rounded.toString()
        else String.format(Locale.ROOT, "%.3f", value).trimEnd('0').trimEnd('.')
    }

    private fun toSignedElevation(value: Double?): Short = when {
        value == null || !value.isFinite() -> NO_DATA.toInt().toShort()
        else -> value.roundToInt().coerceIn(-32767, 32767).toShort()
    }

    private fun toSignedMagnitudeElevation(value: Double?): Short {
        val signed = if (value == null || !value.isFinite()) -32767 else {
            value.roundToInt().coerceIn(-32767, 32767)
        }
        val encoded = if (signed < 0) 0x8000 or -signed else signed
        return encoded.toShort()
    }

    private fun tileName(south: Int, west: Int, extension: String): String =
        "%c%02d%c%03d.%s".format(
            Locale.ROOT,
            if (south < 0) 'S' else 'N',
            kotlin.math.abs(south),
            if (west < 0) 'W' else 'E',
            kotlin.math.abs(west),
            extension,
        )

    private fun dtedHeader(south: Int, west: Int, dimension: Int, intervalTenths: Int): ByteArray {
        val header = ByteArray(DTED_HEADER_SIZE) { ' '.code.toByte() }
        header.putAscii(0, "UHL1")
        header.putAscii(4, dtedAngle(west, longitude = true))
        header.putAscii(12, dtedAngle(south, longitude = false))
        header.putAscii(20, intervalTenths.toString().padStart(4, '0'))
        header.putAscii(24, intervalTenths.toString().padStart(4, '0'))
        header.putAscii(28, "NA  ")
        header.putAscii(32, "U  ")
        header.putAscii(47, dimension.toString().padStart(4, '0'))
        header.putAscii(51, dimension.toString().padStart(4, '0'))
        header.putAscii(80, "DSI")
        header.putAscii(80 + 648, "ACC")
        return header
    }

    private fun dtedAngle(value: Int, longitude: Boolean): String = buildString {
        append(kotlin.math.abs(value).toString().padStart(if (longitude) 3 else 3, '0'))
        append("0000")
        append(if (longitude) if (value < 0) 'W' else 'E' else if (value < 0) 'S' else 'N')
    }

    private fun ByteArray.putAscii(offset: Int, text: String) {
        text.toByteArray(StandardCharsets.US_ASCII).copyInto(this, offset)
    }

    private fun RandomAccessFile.writeIfdEntry(tag: Int, type: Int, count: Long, valueOrOffset: Long) {
        writeShortLe(tag)
        writeShortLe(type)
        writeIntLe(count)
        writeIntLe(valueOrOffset)
    }

    private fun RandomAccessFile.writeShortLe(value: Int) {
        write(value and 0xFF)
        write((value ushr 8) and 0xFF)
    }

    private fun RandomAccessFile.writeIntLe(value: Long) {
        require(value in 0..UINT32_MAX) { "TIFF-Wert außerhalb des UInt32-Bereichs: $value" }
        write((value and 0xFF).toInt())
        write((value ushr 8 and 0xFF).toInt())
        write((value ushr 16 and 0xFF).toInt())
        write((value ushr 24 and 0xFF).toInt())
    }

    private fun RandomAccessFile.writeDoubleLe(value: Double) {
        val bits = java.lang.Double.doubleToRawLongBits(value)
        repeat(8) { shift -> write((bits ushr (shift * 8) and 0xFF).toInt()) }
    }

    private fun align(value: Long, alignment: Long): Long =
        ((value + alignment - 1L) / alignment) * alignment

    private fun Path.partial(): Path = resolveSibling("${fileName}.partial")

    private fun moveAtomically(source: Path, target: Path) {
        try {
            Files.move(source, target, StandardCopyOption.ATOMIC_MOVE)
        } catch (_: AtomicMoveNotSupportedException) {
            Files.move(source, target)
        }
    }

    private fun deleteDirectory(directory: Path) {
        Files.walk(directory).use { paths ->
            paths.sorted(Comparator.reverseOrder()).forEach(Files::deleteIfExists)
        }
    }

    private const val NO_DATA = -32768f
    private const val MAX_GEOTIFF_STRIP_BYTES = 8 * 1024 * 1024
    private const val UINT32_MAX = 0xFFFF_FFFFL
    private const val TIFF_ASCII = 2
    private const val TIFF_SHORT = 3
    private const val TIFF_LONG = 4
    private const val TIFF_DOUBLE = 12
    private const val DTED_HEADER_SIZE = 80 + 648 + 2700
}

private data class OutputGrid(
    val bounds: GeoBounds,
    val width: Int,
    val height: Int,
    val cellSize: Double,
) {
    val cellCount: Long = width.toLong() * height

    fun longitudeAt(column: Int): Double = bounds.west + (column + 0.5) * cellSize
    fun latitudeAt(row: Int): Double = bounds.north - (row + 0.5) * cellSize

    companion object {
        fun create(bounds: GeoBounds, requestedCellSize: Double): OutputGrid {
            require(requestedCellSize.isFinite() && requestedCellSize > 0.0) {
                "Die Rasterauflösung ist ungültig."
            }
            val width = ceil((bounds.east - bounds.west) / requestedCellSize).toLong().coerceAtLeast(1)
            val height = ceil((bounds.north - bounds.south) / requestedCellSize).toLong().coerceAtLeast(1)
            require(width <= Int.MAX_VALUE && height <= Int.MAX_VALUE && width * height <= Int.MAX_VALUE.toLong()) {
                "Das gewählte Ausgaberaster ist zu groß. Bitte wähle einen kleineren Bereich."
            }
            val adjustedBounds = GeoBounds(
                west = bounds.west,
                south = bounds.north - height * requestedCellSize,
                east = bounds.west + width * requestedCellSize,
                north = bounds.north,
            )
            return OutputGrid(adjustedBounds, width.toInt(), height.toInt(), requestedCellSize)
        }
    }
}

private class MergedElevationSampler(
    private val sources: List<ElevationSource>,
    private val areas: List<GeoSelectionRegion>,
) {
    private val sourceIndex = GeographicBucketIndex(sources, ElevationSource::bounds)
    private val areaIndex = GeographicBucketIndex(areas, GeoSelectionRegion::bounds)

    fun sample(longitude: Double, latitude: Double): Double? {
        if (areaIndex.candidates(longitude, latitude).none { area ->
                area.bounds.contains(longitude, latitude) && area.contains(latitude, longitude)
            }
        ) {
            return null
        }
        return sampleSelected(longitude, latitude)
    }

    fun sampleSelected(longitude: Double, latitude: Double): Double? =
        sourceIndex.candidates(longitude, latitude).firstNotNullOfOrNull { source ->
            source.takeIf { it.bounds.contains(longitude, latitude) }?.sample(longitude, latitude)
        }

    fun selectedLongitudeIntervals(latitude: Double): List<ClosedFloatingPointRange<Double>> {
        val intervals = areas.mapNotNull { area ->
            when (area) {
                is GeoRectangleRegion -> area.bounds.takeIf { latitude in it.south..it.north }
                    ?.let { it.west..it.east }
                is GeoCircleRegion -> area.longitudeIntervalAt(latitude)
            }
        }.sortedBy { it.start }
        if (intervals.size < 2) return intervals
        return buildList {
            var start = intervals.first().start
            var end = intervals.first().endInclusive
            intervals.drop(1).forEach { interval ->
                if (interval.start > end) {
                    add(start..end)
                    start = interval.start
                    end = interval.endInclusive
                } else {
                    end = maxOf(end, interval.endInclusive)
                }
            }
            add(start..end)
        }
    }

    private fun GeoCircleRegion.longitudeIntervalAt(latitude: Double): ClosedFloatingPointRange<Double>? {
        if (latitude !in bounds.south..bounds.north) return null
        val centerLatitudeRadians = centerLatitude * PI / 180.0
        val latitudeRadians = latitude * PI / 180.0
        val angularRadius = radiusKilometers / EARTH_MEAN_RADIUS_KM
        val denominator = cos(centerLatitudeRadians) * cos(latitudeRadians)
        if (kotlin.math.abs(denominator) < 0.000_000_000_001) return bounds.west..bounds.east
        val cosine = (cos(angularRadius) - sin(centerLatitudeRadians) * sin(latitudeRadians)) / denominator
        if (cosine > 1.0) return null
        val longitudeDelta = acos(cosine.coerceIn(-1.0, 1.0)) * 180.0 / PI
        return maxOf(bounds.west, centerLongitude - longitudeDelta)..
            minOf(bounds.east, centerLongitude + longitudeDelta)
    }
}

private class GeographicBucketIndex<T>(
    items: List<T>,
    boundsOf: (T) -> GeoBounds,
) {
    private val buckets: Map<Long, List<T>> = buildMap {
        items.forEach { item ->
            val bounds = boundsOf(item)
            for (latitudeBucket in bucket(bounds.south)..bucket(bounds.north)) {
                for (longitudeBucket in bucket(bounds.west)..bucket(bounds.east)) {
                    val key = key(longitudeBucket, latitudeBucket)
                    val bucketItems = get(key)?.toMutableList() ?: mutableListOf()
                    bucketItems += item
                    put(key, bucketItems)
                }
            }
        }
    }

    fun candidates(longitude: Double, latitude: Double): List<T> =
        buckets[key(bucket(longitude), bucket(latitude))].orEmpty()

    private companion object {
        fun bucket(value: Double): Int = floor(value).toInt()

        fun key(longitudeBucket: Int, latitudeBucket: Int): Long =
            (latitudeBucket.toLong() shl 32) xor (longitudeBucket.toLong() and 0xFFFF_FFFFL)
    }
}

internal interface ElevationSource : Closeable {
    val bounds: GeoBounds
    val cellSizeDegrees: Double
    fun sample(longitude: Double, latitude: Double): Double?
}

private fun GeoBounds.contains(longitude: Double, latitude: Double): Boolean =
    longitude in west..east && latitude in south..north

private class HgtElevationSource(
    path: Path,
    override val bounds: GeoBounds,
    private val dimension: Int,
) : ElevationSource {
    private val channel = FileChannel.open(path, StandardOpenOption.READ)
    private val data: MappedByteBuffer = channel.map(FileChannel.MapMode.READ_ONLY, 0, path.fileSize()).apply {
        order(ByteOrder.BIG_ENDIAN)
    }
    override val cellSizeDegrees: Double = 1.0 / (dimension - 1)

    override fun sample(longitude: Double, latitude: Double): Double? {
        val column = ((longitude - bounds.west) / cellSizeDegrees).roundToInt().coerceIn(0, dimension - 1)
        val row = ((bounds.north - latitude) / cellSizeDegrees).roundToInt().coerceIn(0, dimension - 1)
        val value = data.getShort((row * dimension + column) * 2).toInt()
        return value.takeUnless { it == -32768 }?.toDouble()
    }

    override fun close() = channel.close()
}

private class DtedElevationSource(
    path: Path,
    override val bounds: GeoBounds,
    private val width: Int,
    private val height: Int,
) : ElevationSource {
    private val channel = FileChannel.open(path, StandardOpenOption.READ)
    private val data: MappedByteBuffer = channel.map(FileChannel.MapMode.READ_ONLY, 0, path.fileSize()).apply {
        order(ByteOrder.BIG_ENDIAN)
    }
    private val longitudeStep = (bounds.east - bounds.west) / (width - 1)
    private val latitudeStep = (bounds.north - bounds.south) / (height - 1)
    private val recordSize = 8 + height * 2 + 4
    override val cellSizeDegrees: Double = max(longitudeStep, latitudeStep)

    override fun sample(longitude: Double, latitude: Double): Double? {
        val column = ((longitude - bounds.west) / longitudeStep).roundToInt().coerceIn(0, width - 1)
        val row = ((latitude - bounds.south) / latitudeStep).roundToInt().coerceIn(0, height - 1)
        val offset = DTED_DATA_OFFSET + column * recordSize + 8 + row * 2
        require(offset + 2 <= data.limit()) { "Die DTED-Datensätze sind unvollständig." }
        val raw = data.getShort(offset).toInt() and 0xFFFF
        val value = if (raw and 0x8000 != 0) -(raw and 0x7FFF) else raw
        return value.takeUnless { it == -32767 }?.toDouble()
    }

    override fun close() = channel.close()

    private companion object {
        const val DTED_DATA_OFFSET = 80 + 648 + 2700
    }
}

private class AsciiGridElevationSource(path: Path) : ElevationSource {
    private val file = RandomAccessFile(path.toFile(), "r")
    private val width: Int
    private val height: Int
    private val cellSize: Double
    private val noData: Double?
    private val rowOffsets: LongArray
    override val bounds: GeoBounds
    override val cellSizeDegrees: Double get() = cellSize
    private var cachedRow = -1
    private var cachedValues = DoubleArray(0)

    init {
        val values = linkedMapOf<String, String>()
        var dataStart = -1L
        while (true) {
            val position = file.filePointer
            val line = file.readLine() ?: break
            val parts = line.trim().split(Regex("\\s+"), limit = 2)
            val key = parts.firstOrNull()?.lowercase()
            if (parts.size == 2 && key in HEADER_KEYS) {
                values[requireNotNull(key)] = parts[1].trim()
            } else if (line.isNotBlank()) {
                dataStart = position
                break
            }
        }
        width = values.requireInt("ncols")
        height = values.requireInt("nrows")
        cellSize = values.requireDouble("cellsize")
        val west = values["xllcorner"]?.toDoubleOrNull()
            ?: values.requireDouble("xllcenter") - cellSize / 2.0
        val south = values["yllcorner"]?.toDoubleOrNull()
            ?: values.requireDouble("yllcenter") - cellSize / 2.0
        noData = values["nodata_value"]?.toDoubleOrNull()
        bounds = GeoBounds(west, south, west + width * cellSize, south + height * cellSize)
        require(dataStart >= 0L) { "Das ASCII Grid enthält keine Rasterwerte." }
        rowOffsets = LongArray(height)
        file.seek(dataStart)
        for (row in 0 until height) {
            rowOffsets[row] = file.filePointer
            require(file.readLine() != null) { "Das ASCII Grid enthält zu wenige Rasterzeilen." }
        }
    }

    @Synchronized
    override fun sample(longitude: Double, latitude: Double): Double? {
        val column = floor((longitude - bounds.west) / cellSize).toInt().coerceIn(0, width - 1)
        val row = floor((bounds.north - latitude) / cellSize).toInt().coerceIn(0, height - 1)
        if (row != cachedRow) {
            file.seek(rowOffsets[row])
            cachedValues = requireNotNull(file.readLine())
                .trim()
                .split(Regex("\\s+"))
                .map(String::toDouble)
                .toDoubleArray()
            require(cachedValues.size == width) { "Eine ASCII-Grid-Zeile besitzt nicht $width Werte." }
            cachedRow = row
        }
        val value = cachedValues[column]
        return value.takeUnless { !it.isFinite() || (noData != null && it == noData) }
    }

    override fun close() = file.close()

    private fun Map<String, String>.requireInt(key: String): Int =
        get(key)?.toIntOrNull() ?: error("Im ASCII-Grid-Header fehlt '$key'.")

    private fun Map<String, String>.requireDouble(key: String): Double =
        get(key)?.toDoubleOrNull() ?: error("Im ASCII-Grid-Header fehlt '$key'.")

    private companion object {
        val HEADER_KEYS = setOf(
            "ncols", "nrows", "xllcorner", "yllcorner", "xllcenter", "yllcenter", "cellsize", "nodata_value",
        )
    }
}

private class GeoTiffElevationSource(path: Path, override val bounds: GeoBounds) : ElevationSource {
    private val store: DataStore = DataStores.open(path.toFile())
    private val resource: GridCoverageResource = (store as? Aggregate)
        ?.components()
        ?.filterIsInstance<GridCoverageResource>()
        ?.firstOrNull()
        ?: error("Die GeoTIFF-Datei enthält kein Höhenraster.")
    private val coverage = resource.read(null).forConvertedValues(true)
    private val geometry = coverage.gridGeometry
    private val extent = geometry.extent
    private val image: RenderedImage = coverage.render(null)
    private val sourceCrs = geometry.coordinateReferenceSystem
        ?: error("Das GeoTIFF besitzt kein Koordinatenreferenzsystem.")
    private val wgs84ToSource: MathTransform = CRS.findOperation(
        CommonCRS.WGS84.normalizedGeographic(),
        sourceCrs,
        null,
    ).mathTransform
    private val sourceToGrid: MathTransform = geometry.getGridToCRS(PixelInCell.CELL_CENTER).inverse()
    private val noData = coverage.sampleDimensions.firstOrNull()?.background?.orElse(null)?.toDouble()
        ?: readGeoTiffGdalNoData(path)
    override val cellSizeDegrees: Double = max(
        (bounds.east - bounds.west) / extent.getSize(0),
        (bounds.north - bounds.south) / extent.getSize(1),
    )

    @Synchronized
    override fun sample(longitude: Double, latitude: Double): Double? {
        val geographic = DirectPosition2D(longitude, latitude)
        val source = wgs84ToSource.transform(geographic, null)
        val grid = sourceToGrid.transform(source, null)
        val column = grid.getOrdinate(0).roundToInt().toLong()
        val row = grid.getOrdinate(1).roundToInt().toLong()
        if (!extent.contains(column, row)) return null
        val x = image.minX + (column - extent.getLow(0)).toInt()
        val y = image.minY + (row - extent.getLow(1)).toInt()
        val tileX = floor((x - image.tileGridXOffset).toDouble() / image.tileWidth).toInt()
        val tileY = floor((y - image.tileGridYOffset).toDouble() / image.tileHeight).toInt()
        val value = image.getTile(tileX, tileY).getSampleDouble(x, y, 0)
        return value.takeUnless { !it.isFinite() || (noData != null && it == noData) }
    }

    override fun close() = store.close()
}
