GDAL（一）：插值和裁剪（python）

在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)