GeoTools插值无效：生成GeoTIFF前后无变化问题求助

我有包含lat、lon和z值的空间数据，以及对应的多边形范围。需要创建可通过GeoServer发布的GeoTIFF，要求多边形内缺失值做插值处理实现平滑效果，多边形外设为NoData。但用GeoTools实现时，插值前后的GeoTIFF完全一致，附上Kotlin代码及效果对比图，求排查问题原因。

class GeoToolsExperiment {
    fun execute(pointWithDatas: List<PointWithData>, polygon: Polygon) {

        val crs: CoordinateReferenceSystem = CRS.decode("EPSG:3857")

        val bounds = ReferencedEnvelope(polygon.envelopeInternal, crs)
        val factory: GridCoverageFactory = CoverageFactoryFinder.getGridCoverageFactory(null)

        val width = (bounds.width * 100000).toInt()
        val height = (bounds.height * 100000).toInt()

        val writableRaster: WritableRaster = createRaster(width, height)
        val gc: GridCoverage2D = createGridCoverage(factory, writableRaster, bounds)

        fillRaster(pointWithDatas, gc, writableRaster)
        
        val ci = interpolate(gc)
      
        writeToTif(ci)

    }

    private fun writeToTif(ci: GridCoverage?) {
        val outFile = "test.tif"
        val out = File(outFile)
        val format: AbstractGridFormat = GeoTiffFormat()
        val writer = format.getWriter(out)
        try {
            writer.write(ci)
            writer.dispose()
        } catch (e: IllegalArgumentException) {
            e.printStackTrace()
        } catch (e: IOException) {
            e.printStackTrace()
        }
    }

    private fun interpolate(gc: GridCoverage2D): GridCoverage {
        val interp = Interpolation.getInstance(Interpolation.INTERP_BILINEAR)
        val interpolate = Operations.DEFAULT.interpolate(gc, interp)
        return interpolate as GridCoverage2D
    }

    private fun createGridCoverage(
        factory: GridCoverageFactory,
        writableRaster: WritableRaster,
        bounds: ReferencedEnvelope
    ): GridCoverage2D {
        val gc: GridCoverage2D = factory.create("band1", writableRaster, bounds)
        return gc
    }

    private fun createRaster(width: Int, height: Int): WritableRaster {
        val writableRaster: WritableRaster = RasterFactory.createBandedRaster(
            DataBuffer.TYPE_DOUBLE, width,
            height, 1, null
        )
        return writableRaster
    }

    private fun fillRaster(
        pointWithDatas: List<PointWithData>,
        gc: GridCoverage2D,
        writableRaster: WritableRaster
    ) {
        pointWithDatas.forEach { pwd ->
            val directPosition = DirectPosition2D(pwd.point.x, pwd.point.y)
            val transformPointToGridCoverage = gc.gridGeometry.worldToGrid(directPosition)
            val data = DoubleArray(1)
            data[0] = pwd.data
            writableRaster.setDataElements(transformPointToGridCoverage.x, transformPointToGridCoverage.y, data)
        }
    }
}

多边形范围

矢量数据

代码运行结果

期望效果（QGIS手动生成）

问题核心原因 #

1. 初始栅格无NoData标记
   你创建的栅格默认值为0.0（Double类型默认值），GeoTools的插值操作仅针对NoData区域生效。当前所有栅格单元都被判定为有效数据，插值逻辑不会触发。
2. 插值方式选错
   INTERP_BILINEAR是栅格内部的平滑插值，只能对已有栅格单元的边缘做平滑，无法基于离散点向空白区域扩散数值，和QGIS中基于点的空间插值（如IDW、克里金）不是同一类操作。
3. 未处理多边形范围
   代码未区分多边形内外区域，即使插值生效，也会填充整个外接矩形范围，不符合需求。

修正方案 #

1. 初始化栅格为NoData

修改createRaster方法，将所有单元初始化为Double.NaN（GeoTools默认识别的Double类型NoData值）：

private fun createRaster(width: Int, height: Int): WritableRaster {
    val writableRaster = RasterFactory.createBandedRaster(
        DataBuffer.TYPE_DOUBLE, width, height, 1, null
    )
    // 批量设置所有单元为NaN
    val noData = DoubleArray(1) { Double.NaN }
    for (x in 0 until width) {
        for (y in 0 until height) {
            writableRaster.setDataElements(x, y, noData)
        }
    }
    return writableRaster
}

2. 替换为基于离散点的IDW插值

用反距离加权插值（IDW）替代原有栅格平滑插值，实现离散点的区域填充：

private fun interpolate(gc: GridCoverage2D, polygon: Polygon): GridCoverage2D {
    // 提取有效样本点
    val samples = mutableListOf<Point2D>()
    val values = mutableListOf<Double>()
    val raster = gc.renderedImage.data as Raster
    val gridToWorld = gc.gridGeometry.gridToWorld(PixelInCell.CELL_CENTER)

    for (x in 0 until raster.width) {
        for (y in 0 until raster.height) {
            val value = raster.sampleDouble(x, y, 0)
            if (!value.isNaN()) {
                val worldPos = gridToWorld.transform(DirectPosition2D(x.toDouble(), y.toDouble()), null)
                samples.add(Point2D.Double(worldPos.x, worldPos.y))
                values.add(value)
            }
        }
    }

    // 创建IDW插值器
    val idw = Interpolator2D(Interpolator2D.INVERSE_DISTANCE, samples.toTypedArray(), values.toDoubleArray())
    val newRaster = createRaster(raster.width, raster.height)
    val worldToGrid = gc.gridGeometry.worldToGrid(PixelInCell.CELL_CENTER)

    // 填充插值结果，同时处理多边形范围
    for (x in 0 until newRaster.width) {
        for (y in 0 until newRaster.height) {
            val worldPos = gridToWorld.transform(DirectPosition2D(x.toDouble(), y.toDouble()), null)
            val point = Point(worldPos.x, worldPos.y)
            if (polygon.contains(point)) {
                val interpolatedValue = idw.interpolate(worldPos.x, worldPos.y)
                newRaster.setSample(x, y, 0, interpolatedValue)
            }
            // 多边形外保持NoData（已初始化）
        }
    }

    return CoverageFactoryFinder.getGridCoverageFactory(null)
        .create("interpolated", newRaster, gc.envelope)
}

3. 修改execute方法传递多边形参数

调用插值方法时传入polygon：

val ci = interpolate(gc, polygon)

内容的提问来源于stack exchange，提问作者glennvdv