處理網格（Raster）資料

在 Python 中處理地理空間資料

Joris Van den Bossche

Open source software developer and teacher, GeoPandas maintainer

數值網格示意

圖片來源：QGIS 說明文件

Raster 資料

降雨機率地圖

多波段的 Raster 資料

rasterio 套件

import rasterio

對 GDAL 的「Pythonic」繫結
讀取與寫入 raster 檔案
處理工具（遮罩、重投影、重取樣等）

https://rasterio.readthedocs.io/en/latest/

開啟 raster 檔案

import rasterio

src = rasterio.open("DEM_world.tif")

中繼資料：

src.count

src.width, src.height

(4320, 2160)

Raster 資料 = numpy 陣列

array = src.read()

標準 numpy 陣列：

array

array([[[-4290, -4290, -4290, ..., -4290, -4290, -4290],
        [-4278, -4278, -4278, ..., -4278, -4278, -4278],
        [-4269, -4269, -4269, ..., -4269, -4269, -4269],
        ...,
        [ 2804,  2804,  2804, ...,  2804,  2804,  2804],
        [ 2804,  2804,  2804, ...,  2804,  2804,  2804],
        [ 2804,  2804,  2804, ...,  2804,  2804,  2804]]], dtype=int16)

繪製 raster 資料集

使用 rasterio.plot.show() 方法：

import rasterio.plot

rasterio.plot.show(src, cmap='terrain')

依向量資料擷取資訊

rasterstats：依向量幾何計算地理空間 raster 資料集的摘要統計（https://github.com/perrygeo/python-rasterstats）

用 rasterstats 擷取 raster 數值

點狀向量：

rasterstats.point_query(geometries, "path/to/raster", 
                        interpolation='nearest'|'bilinear')

多邊形向量：

rasterstats.zonal_stats(geometries, "path/to/raster",
                        stats=['min', 'mean', 'max'])

用 rasterstats 擷取 raster 數值

result = rasterstats.zonal_stats(countries.geometry, "DEM_gworld.tif", 
                                 stats=['mean'])

countries['mean_elevation'] = pd.DataFrame(result)

countries.sort_values('mean_elevation', ascending=False).head()

            name    continent                     geometry  mean_elevation
157   Tajikistan         Asia  POLYGON ((74.98 37.41, ...      3103.231105
85    Kyrgyzstan         Asia  POLYGON ((80.25 42.34, ...      2867.717142
24        Bhutan         Asia  POLYGON ((91.69 27.77, ...      2573.559846
119        Nepal         Asia  POLYGON ((81.11 30.18, ...      2408.907816
6     Antarctica   Antarctica  (POLYGON ((-59.57 -80.04...     2374.075028
..           ...          ...                          ...             ...

一起來練習吧！

在 Python 中處理地理空間資料