在GIS中，經常遇到原始資料是點，但是展示的時候點展示並不好，能區域內連續展示最好了。

這個時候就需要用到插值，把點轉成連續的面展示。

大部分的GIS軟體中都有插值的工具可以直接使用，不過當遇到批量插值的時候，工具用起來就比較費時間了。

所以就想到寫程式碼，進行批量操作，這樣可以執行程式碼，就不用管了。

既然用程式碼批量操作，GDAL就是最佳選擇。本想是用 .net 的，因為是引用C++的dll，用起來很不方便，而且文件、資料比較少，就選擇了 python來寫。

一、網格化、演算法簡介

網上找對應的實現程式碼是比較好找，但是對這個實現的基本介紹、引數怎麼設定等基本沒有。

如果只是使用是沒問題，但是想用好還是有點問題的。

這裡說下自己歷程：

• 查詢程式碼實現，這個好找

• 理解引數，這個相對較少

• 演算法配置，網上基本沒有

最後發現是在官網找到的介紹，比較簡單，配置使用是沒問題的。

重點：對於成熟的開源專案，看官網！看官網！看官網！這樣可以避免很多坑。

這裡對於具體的演算法也就不拿過來了，大家直接去官網看：

Grid 簡介 （對應中文版）

Grid 引數解釋、配置（對應中文版）

有一個對插值演算法解釋挺好的，在這裡貼出來：（原文地址）

反距離權重插值

'invdist:power=3.6:smoothing=0.2:radius1=0.0:radius2=0.0:angle=0.0:max_points=0:min_points=0:nodata=0.0
 
'

引數解釋：

• invdist：表示演算法名稱 - 反距離權重，即距離越近權重越大，對待計算的網格影像越大；

• power：距離對權重的影響程度）（數字表示指數），預設值為2, 值越大表示距離越近的點對插值的影響程指數倍增；對於該引數的設定，若點較密集且均勻分佈，則值設定在2以下，若點相對較少且分佈不均，則為了保障插值效果，可將引數設定在3以上；

• smoothing：平滑係數，預設為0，數值越大表示越平滑，同時精度也會越低；

• radius1：表示橢圓x軸方向上的半徑；

• radius2：表示橢圓y軸方向上的半徑；

• angle：表示橢圓旋轉的弧度；

• max_points：表示參與插值的最大點數，預設值0表示搜尋到多少就是多少；

• min_points：表示參與插值的最小點數，預設值0表示搜尋到多少就是多少；

• nodata：對nodata的填充值

二、程式碼實現

程式碼實現起來還是挺簡單的：

# field 這裡是要插值的欄位
def insertRaster(field):
    startTime = ps.to_datetime(datetime.datetime.now())
    print('開始插值：%s' % startTime)
    point_file = 'xxx.shp'
    output_file = 'xxx' + field + '.tif'

    opts = gdal.GridOptions(format="GTiff", outputType=gdal.GDT_Float32, width=2472, height=1000,
                            algorithm="invdist:power=3.5:smothing=0.0:radius=1.0:max_points=12:min_points=0:nodata=0.0", zfield=field)

    gdal.Grid(destName=output_file, srcDS=point_file, options=opts)
    endTime = ps.to_datetime(datetime.datetime.now())
    useTime = endTime - startTime
    print('插值完成！用時：%s s' % useTime.total_seconds())

這裡面的 GridOptions 裡面的引數都可以看文件，演算法是上面給的連結，可以根據裡面進行配置。

三、裁剪

為什麼要在這裡帶上裁剪。

插值的輸出結果，都是一個矩形。大部分場景下，這和我們需要使用的範圍不一致。而且源資料的範圍也不一定是矩形，在源資料範圍以外，插值的結果都是無效的，所以要裁減掉。

裁剪對應的簡介、引數配置也都可以在官網找到。對應直接上程式碼：

def cutRaster():
    input_shp = 'xxx.shp'
    input_raster = 'xxx-r '.tif'
    result_raster = 'xxx-r'_c.tiff'

    if os.path.exists(input_raster):
        startTime = ps.to_datetime(datetime.datetime.now())
        print('開始裁剪：%s' % startTime)
        ds = gdal.Warp(result_raster, input_raster, format='GTiff',
                       cutlineDSName=input_shp, dstNodata=0)
        endTime = ps.to_datetime(datetime.datetime.now())
        useTime = endTime - startTime
        print('裁剪完成，用時：%s s' % useTime.total_seconds())
    else:
        print('檔案不存在：%s' % input_raster)