用Python進行柵格數據的分區統計和批量提取
有時候我們會有這樣的想法,就是針對某個區域的柵格數據,要提取它的平均值或者其他統計指標,比如在一個省內提取多年的降雨數據,最后分區域地計算一些統計值,或者從多個柵格數據中提取某個區域的數值形成一個序列。為了方便,畫一個示意圖看看,比如就像提取這個區域中的某一個市的區域,然后形成一個序列數據,這就可以使用rasterstats庫了,此外的分區統計也可以用這個庫
這個實驗使用的數據格式分別是柵格(*.tif)和矢量(.shp),之后的分區統計操作和柵格數據的提取都是源于這兩類數據。為了能使用上這個rasterstats庫,選擇了在google colab平臺運行腳本,因為安裝庫實在是太方便了,在win上老是安裝不上的,在google notebook立馬就搞定了,而且可以把數據存儲到谷歌云盤,直接在notebook中就是可以鏈接使用的
那么現在就開始做測試,使用的數據就是左側的柵格和矢量數據集導入相關的模塊
import geopandas as gpdimport pandas as pdimport numpy as npimport matplotlib.pyplot as pltimport rasterioimport rasterstatsfrom rasterio.plot import show# show()方法用來展示柵格圖形from rasterio.plot import show_hist# 用來展示直方圖import cartopy.crs as ccrsimport cartopy.feature as cfeaturefrom cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
使用geopandas和rasterio分別讀取矢量和柵格數據
# 使用geopandas讀取矢量數據districts = gpd.read_file(’/content/drive/MyDrive/Datashpraster/Data/Districts/districts.shp’)# 使用rasterio讀取柵格數據,柵格數據和矢量數據的坐標投影需要一致raster = rasterio.open(’/content/drive/MyDrive/Datashpraster/Data/Rainfall Data Rasters/2020-4-1.tif’)
# 把矢量數據和柵格數據繪制到一個axis上,這個axis不是坐標軸,而是圖形plt.rcParams[’font.family’] = ’Times New Roman’plt.rcParams[’font.size’] = 20fig, (ax1,ax2) = plt.subplots(1,2,figsize=(15,6))show(raster, ax=ax1,title=’Rainfall’)# 讀取進來的矢量數據可以直接調用gpd的plot()方法繪制districts.plot(ax=ax1, facecolor=’None’, edgecolor=’red’)show_hist(raster,ax=ax2,title=’hist’)plt.show()
先繪制一下結果看看
讀取柵格數據:
# 提取雨量柵格值到numpy數組# 遵循GDAL規則從第一波段讀取rainfall_data = raster.read(1)rainfall_data
開始分區統計:
# 設置坐標變換信息affine = raster.transform# 準備開始進行空間分區計算# 第一個參數是矢量分區,第二個是柵格,第三個是坐標變換信息,第四個是統計均值avg_rallrain = rasterstats.zonal_stats(districts,rainfall_data,affine=affine,stats=[’mean’],geojson_out=True)# avg_rallrain# 除了統計平均值之外,還有最大最小值那些
繪制一下,只是一個簡單的圖形而已
當然第二部分更有意思,就是從多個分散的柵格數據中提取數據形成一個序列
,就是這些tif數據
loop這些柵格數據集:
獲得提取到的結果,沒錯,就是這么一個序列數據,然后就是繪圖了
轉換數據格式
# 將Date列轉為時間型data[’Date’] = pd.to_datetime(data[’Date’], infer_datetime_format=True)# print(data)data[’Date’] = data[’Date’].dt.dateprint(data)
繪圖結果就是簡單的圖形而已
# 準備繪制圖形fig,(ax1,ax2)= plt.subplots(2,1,figsize=(18,6))plt.rcParams[’font.size’] = 15data.plot(x=’Date’, y=’Average_RF_Porto’, ax=ax1, kind=’bar’, title=’Avg_Rail_Porto’)data.plot(x=’Date’, y=’Average_RF_Faro’, ax=ax2, kind=’bar’, title=’Avg_Rail_Faro’,color=’red’)#自動調整圖形的分布plt.tight_layout()plt.show()
結果就這樣一個序列圖,目的就是從柵格提取指定的研究區,然后提取柵格的值,再來繪圖
雖然感覺不是那么花里胡哨的圖,但這個應該還是比較實用的,特別是大批量提取柵格值的時候。由于在google colab里面操作的步驟比較多,中間可能有省略的地方,但重要的應該都在文中了,當然也可以遷移運用到其他地方,也可以查看一下這個第三方庫的教程,比如read(1)是什么意思,官網的docs就寫得有,實在是很方便的
以上就是用Python進行柵格數據的分區統計和批量提取的詳細內容,更多關于Python 柵格數據的分區統計和批量提取 的資料請關注好吧啦網其它相關文章!
相關文章: