GEE OEEL kriging 包Python API版本（geemap）-克里金插值

背景

众所周知OEEL包是一个好用的GEE libaray, 根据吴秋生老师和OEEL官网的教程，均表示使用下列代码，可以直接在python api中调用JavaScript写的OEEL函数。

oeel = geemap.requireJS()

然而我试了好久好久，一直都是在加载的状态....

所以我就把用到的OEEL函数克里金补缺失值的函数源码转为了Python版本。

代码

#源码链接：https://github.com/open-geocomputing/OpenEarthEngineLibrary/blob/master/Image/kriging

#Python版本
# SPDX-License-Identifier: LGPL-3.0-or-later

import geemap
import ee

# 定义输入参数
inputs = [
    {'name': 'covFun', 'description': 'The covariance function, fun(ee.Array)', 'type': 'function', 'optional': False},
    {'name': 'radius', 'description': 'The radius of the windows used to get neighbors.', 'type': 'integer', 'default_value': 10, 'optional': False},
    {'name': 'skipInfomed', 'description': 'Do not compute on informed locations. To set to false in presence of a kernel with a 0 in the center', 'type': 'boolean', 'default_value': True, 'optional': False},
    {'name': 'kernel', 'description': 'Neighbors with weight equal to 0 are ignored, the size needs to be compatible with the radius.', 'type': 'ee.Kernel|ee.Array', 'optional': True},
    {'name': 'image', 'description': 'The image on which to run the computation.', 'type': 'ee.Image', 'optional': True},
    {'name': 'imageCollection', 'description': 'The image collection on which to run the computation', 'type': 'ee.ImageCollection', 'optional': True},
    {'name': 'Return', 'description': 'Return', 'type': 'function|ee.Image', 'default_value': None, 'optional': True}
]

reference = {
    'name': 'kriging',
    'license': 'LGPL-3.0-or-later',
    'description': 'Compute a gap filling interpolation using kriging. It returns a function, or apply on an image if "image" is informed and return the result of the computation.',
    'contributors': ['Mathieu Gravey']
}

def create_documentation():
    return {'inputs': inputs, 'reference': reference}

def create_function(covFun, radius=10, skipInfomed=True, kernel=None, image=None, imageCollection=None):
    halfSize = radius
    dataKernel = kernel

    X = ee.Array(ee.List.sequence(-halfSize, halfSize, 1)).repeat(1, halfSize * 2 + 1)
    Y = ee.Array(ee.List.sequence(-halfSize, halfSize, 1)).repeat(1, halfSize * 2 + 1).transpose()
    
    numberOfElement = (halfSize * 2 + 1) * (halfSize * 2 + 1)
    
    dist = ee.Array.cat([X, Y], 2)
    
    distK = dist.reshape([-1, 1, 2]).repeat(1, numberOfElement).subtract(dist.reshape([1, -1, 2]).repeat(0, numberOfElement))
    distK = distK.pow(2).reduce(ee.Reducer.sum(), [2]).sqrt().reshape([(halfSize * 2 + 1) * (halfSize * 2 + 1), (halfSize * 2 + 1) * (halfSize * 2 + 1)])
    K = covFun(distK)
    K = K.add(ee.Array.identity(numberOfElement)).pad([(numberOfElement + 1), (numberOfElement + 1)], 1).subtract(ee.Array.identity((numberOfElement + 1)))
    
    K0 = covFun(dist.pow(2).reduce(ee.Reducer.sum(), [2]).sqrt().reshape([-1, 1])).pad([numberOfElement + 1, 0], 1)
    
    ImageKc = ee.Image.constant(K)
    ImageK0c = ee.Image.constant(K0)
    
    kernel = ee.Kernel.square(halfSize)
    
    if dataKernel and isinstance(dataKernel, ee.Kernel):
        ker = dataKernel.getInfo()['weights'].replace('\n', ',')
        val = '[' + ker[1:-1] + ']'
        dataKernel = ee.Array(json.loads(val))
        print(dataKernel)
    def fun(im):
        mask = im.mask().neighborhoodToArray(kernel=kernel, defaultValue=0)
        if dataKernel:
            mask = mask.multiply(dataKernel)
            
        mask = mask.arrayReshape(ee.Array([-1]), 1).arrayPad([numberOfElement + 1], 1)
        r = im.unmask(0).neighborhoodToArray(kernel=kernel, defaultValue=0).arrayReshape(ee.Array([-1]), 1).arrayPad([numberOfElement + 1], 0).arrayMask(mask)
        
        mask = mask.arrayReshape(ee.Array([-1, 1]), 2)
        
        ImageK0 = ImageK0c.arrayMask(mask)
        
        dataSize = ImageK0.gt(0).arrayReduce(ee.Reducer.sum(), [0]).arrayGet([0, 0]).gt(3)
        
        coef = ImageKc.arrayMask(mask).arrayTranspose().arrayMask(mask).matrixSolve(ImageK0).arrayReshape(ee.Array([-1]), 1)
        estimate = coef.multiply(r).arrayReduce(ee.Reducer.sum(), [0]).arrayGet([0]).updateMask(dataSize)
        variance = coef.multiply(ImageK0.arrayReshape(ee.Array([-1]), 1)).arrayReduce(ee.Reducer.sum(), [0]).subtract(covFun(ee.Array([0]))).arrayGet([0]).subtract(coef.arrayGet([-1]).multiply(2)).multiply(-1).updateMask(dataSize)
        
        if skipInfomed:
            estimate = estimate.where(im, im)
            variance = variance.where(im, 0)
        
        result = im.addBands([estimate.rename('estimate'), variance.rename('variance')])
        
        result = result.reproject(im.projection())
        return result
    
    if image:
        return fun(image)
    if imageCollection:
        return imageCollection.map(fun)
    
    return fun

#调用
map = geemap.Map()
someData = ee.Image("NASA/JPL/global_forest_canopy_height_2005")
def covFun(dist):
    return dist.multiply(-0.1).exp().multiply(140)
result = create_function(covFun=covFun, radius=10, image=someData)
map.addLayer(someData, {'palette': ['green']}, 'someData')
map.addLayer(result.select('estimate'), {'palette': ['red','green','blue']}, 'Kriging Result')
geometry = ee.Geometry.Point([10.44104168672001, 45.596877548727065])
map.centerObject(geometry,10)
map