利用Pyhton对数据分区统计

分区统计是一个用于计算另一个数据集的区域内栅格数据值的统计数据，包括均值，中位数，最大值，最小值，总和等。本章内容将从TIF数据中批量提取SHP区域范围统计值。

涉及的库包括glob、geopandas、rasterio、rasterstats

导入库

import glob
import numpy as np
import rasterio
import rasterstats
import pandas as pd

读取TIF文件目录

[继续安利，glob是一个很好用的读取文件夹中特定文件的库]

Input_folder = r'.\tif\output'
# 读取所有nc数据
data_list = glob.glob(Input_folder + '/*.tif')
data_list[:5]

分区统计

mean = []
for i in data_list:
    mean.append(rasterstats.zonal_stats('shp/Uzb.shp', i))

out[3]:
[[{'min': 1.6692631244659424,
   'max': 32.75654220581055,
   'mean': 7.104131728566754,
   'count': 764}],
 [{'min': 6.693352699279785,
   'max': 52.895751953125,
   'mean': 17.020224046956805,
   'count': 764}],
 [{'min': 15.254048347473145,
   'max': 61.91945266723633,
   'mean': 29.42128701734293,
   'count': 764}],
 [{'min': 9.328380584716797,
   'max': 94.33841705322266,
   'mean': 31.399040821335078,
   'count': 764}],
 [{'min': 3.756659984588623,
   'max': 117.5350112915039,
   'mean': 25.642033601930628,
   'count': 764}],
 [{'min': 1.7687196731567383,
   'max': 112.15010833740234,
   'mean': 16.417166025850786,
   'count': 764}],
 [{'min': 1.5461831092834473,
   'max': 104.57400512695312,
   'mean': 14.721388252617801,
   'count': 764}],
 [{'min': 1.0155596733093262,
   'max': 65.43617248535156,
   'mean': 8.65694854896106,
   'count': 764}],
 [{'min': 0.7461696267127991,
   'max': 20.19913101196289,
   'mean': 5.3231859456806285,
   'count': 764}],
 [{'min': 0.38746243715286255,
   'max': 33.00067138671875,
   'mean': 6.337624754581152,
   'count': 764}],
 [{'min': 0.357392281293869,
   'max': 23.879793167114258,
   'mean': 6.345281311354712,
   'count': 764}],
 [{'min': 1.8076934814453125,
   'max': 26.18817138671875,
   'mean': 8.100336044870746,
   'count': 764}]]

默认的策略是包括沿直线渲染路径的所有像素（对于直线），或者中心点在多边形内的单元（对于多边形）。另外也可以选择all_touched策略，通过包括它所接触的所有像素来栅格化几何体。

左边是默认False，右边是True状态。

onal_stats('shp/Uzb.shp', i, all_touched=True)

根据官方文档，也可以设置需要的统计函数。

def mymean(x): 
    return np.ma.mean(x)

mean = []
for i in data_list:
    mean.append(rasterstats.zonal_stats('shp/Uzb.shp', i, stats="count", add_stats={'mymean':mymean}))
    
mean

out[4]:
[[{'count': 764, 'mymean': 7.104131728566754}],
 [{'count': 764, 'mymean': 17.020224046956805}],
 [{'count': 764, 'mymean': 29.42128701734293}],
 [{'count': 764, 'mymean': 31.399040821335078}],
 [{'count': 764, 'mymean': 25.642033601930628}],
 [{'count': 764, 'mymean': 16.417166025850786}],
 [{'count': 764, 'mymean': 14.721388252617801}],
 [{'count': 764, 'mymean': 8.65694854896106}],
 [{'count': 764, 'mymean': 5.3231859456806285}],
 [{'count': 764, 'mymean': 6.337624754581152}],
 [{'count': 764, 'mymean': 6.345281311354712}],
 [{'count': 764, 'mymean': 8.100336044870746}]]

sum(mymean)
out[5]:
176.4886481000491

与arcgis的区域统计结果一致