利用Pyhton对气象数据分区统计

josephsun

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

1. import glob
   
2. import numpy as np
   
3. import rasterio
   
4. import rasterstats
   
5. import pandas as pd

复制代码

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

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

复制代码

分区统计

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

复制代码

1. out[3]:
   
2. [[{'min': 1.6692631244659424,
   
3.    'max': 32.75654220581055,
   
4.    'mean': 7.104131728566754,
   
5.    'count': 764}],
   
6. [{'min': 6.693352699279785,
   
7.    'max': 52.895751953125,
   
8.    'mean': 17.020224046956805,
   
9.    'count': 764}],
   
10. [{'min': 15.254048347473145,
    
11.    'max': 61.91945266723633,
    
12.    'mean': 29.42128701734293,
    
13.    'count': 764}],
    
14. [{'min': 9.328380584716797,
    
15.    'max': 94.33841705322266,
    
16.    'mean': 31.399040821335078,
    
17.    'count': 764}],
    
18. [{'min': 3.756659984588623,
    
19.    'max': 117.5350112915039,
    
20.    'mean': 25.642033601930628,
    
21.    'count': 764}],
    
22. [{'min': 1.7687196731567383,
    
23.    'max': 112.15010833740234,
    
24.    'mean': 16.417166025850786,
    
25.    'count': 764}],
    
26. [{'min': 1.5461831092834473,
    
27.    'max': 104.57400512695312,
    
28.    'mean': 14.721388252617801,
    
29.    'count': 764}],
    
30. [{'min': 1.0155596733093262,
    
31.    'max': 65.43617248535156,
    
32.    'mean': 8.65694854896106,
    
33.    'count': 764}],
    
34. [{'min': 0.7461696267127991,
    
35.    'max': 20.19913101196289,
    
36.    'mean': 5.3231859456806285,
    
37.    'count': 764}],
    
38. [{'min': 0.38746243715286255,
    
39.    'max': 33.00067138671875,
    
40.    'mean': 6.337624754581152,
    
41.    'count': 764}],
    
42. [{'min': 0.357392281293869,
    
43.    'max': 23.879793167114258,
    
44.    'mean': 6.345281311354712,
    
45.    'count': 764}],
    
46. [{'min': 1.8076934814453125,
    
47.    'max': 26.18817138671875,
    
48.    'mean': 8.100336044870746,
    
49.    'count': 764}]]

复制代码

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

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

复制代码

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

1. def mymean(x):
   
2.     return np.ma.mean(x)
   
3. 
   
4. mean = []
   
5. for i in data_list:
   
6.     mean.append(rasterstats.zonal_stats('shp/Uzb.shp', i, stats="count", add_stats={'mymean':mymean}))
   
7.    
   
8. mean

复制代码

1. out[4]:
   
2. [[{'count': 764, 'mymean': 7.104131728566754}],
   
3. [{'count': 764, 'mymean': 17.020224046956805}],
   
4. [{'count': 764, 'mymean': 29.42128701734293}],
   
5. [{'count': 764, 'mymean': 31.399040821335078}],
   
6. [{'count': 764, 'mymean': 25.642033601930628}],
   
7. [{'count': 764, 'mymean': 16.417166025850786}],
   
8. [{'count': 764, 'mymean': 14.721388252617801}],
   
9. [{'count': 764, 'mymean': 8.65694854896106}],
   
10. [{'count': 764, 'mymean': 5.3231859456806285}],
    
11. [{'count': 764, 'mymean': 6.337624754581152}],
    
12. [{'count': 764, 'mymean': 6.345281311354712}],
    
13. [{'count': 764, 'mymean': 8.100336044870746}]]

复制代码

1. sum(mymean)
   
2. out[5]:
   
3. 176.4886481000491

复制代码

与arcgis的区域统计结果一致

VV09

不错不错，学习了~~

sms

学习学习一下