Python中怎么对栅格数据做提取和统计操作
Admin 2022-07-29 群英技术资讯 1675 次浏览
这篇文章给大家介绍了“Python中怎么对栅格数据做提取和统计操作”的相关知识,讲解详细,步骤过程清晰,有一定的借鉴学习价值,因此分享给大家做个参考,感兴趣的朋友接下来一起跟随小编看看吧。有时候我们会有这样的想法,就是针对某个区域的栅格数据,要提取它的平均值或者其他统计指标,比如在一个省内提取多年的降雨数据,最后分区域地计算一些统计值,或者从多个栅格数据中提取某个区域的数值形成一个序列。为了方便,画一个示意图看看,比如就像提取这个区域中的某一个市的区域,然后形成一个序列数据,这就可以使用rasterstats库了,此外的分区统计也可以用这个库
这个实验使用的数据格式分别是栅格(*.tif)和矢量(.shp),之后的分区统计操作和栅格数据的提取都是源于这两类数据。为了能使用上这个rasterstats库,选择了在google colab平台运行脚本,因为安装库实在是太方便了,在win上老是安装不上的,在google notebook立马就搞定了,而且可以把数据存储到谷歌云盘,直接在notebook中就是可以链接使用的

那么现在就开始做测试,使用的数据就是左侧的栅格和矢量数据集
导入相关的模块
import geopandas as gpd import pandas as pd import numpy as np import matplotlib.pyplot as plt import rasterio import rasterstats from rasterio.plot import show # show()方法用来展示栅格图形 from rasterio.plot import show_hist # 用来展示直方图 import cartopy.crs as ccrs import cartopy.feature as cfeature from 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'] = 20 fig, (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.date print(data)

绘图结果就是简单的图形而已
# 准备绘制图形 fig,(ax1,ax2)= plt.subplots(2,1,figsize=(18,6)) plt.rcParams['font.size'] = 15 data.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就写得有,实在是很方便的
免责声明:本站发布的内容(图片、视频和文字)以原创、转载和分享为主,文章观点不代表本网站立场,如果涉及侵权请联系站长邮箱:mmqy2019@163.com进行举报,并提供相关证据,查实之后,将立刻删除涉嫌侵权内容。
猜你喜欢
python 中使用open打开某个文件写入时,往往会发现需要写入的文件不在同级目录下。这样就需要根据文件的路径来找到并打开。但往往有时绝对路径和相对路径,写入不正确就会打开失败。
这篇文章主要介绍了Python assert断言声明,遇到错误则立即返回问题,具有很好的参考价值,希望对大家有所帮助。如有错误或未考虑完全的地方,望不吝赐教
弹出框是自动化测试中一种常见的元素,这种元素通常是客户端自带的,下面这篇文章主要给大家介绍了关于Selenium定位浏览器弹窗方法的相关资料,文中通过实例代码介绍的非常详细,需要的朋友可以参考下
常规版本简单的for循环遍历x_n=["x1","x2","x3"]forxinx_n:print(x)>>x1x2x3借助range()和len()x_n=["x1","x2","x3"]foriinrange(len(x_n)):print(...
property是一个类,可以把一个方法当做属性进行使用,这样做可以简化代码使用。实际上就是装饰类中属性的gettersetter方法,使得属性可以通过对象.属性的方式获取或设置 使用property的两种方式装饰器方式类属性方式2.装饰器方式@property修饰获取的方法getter,方法名必须和属性名一样@age.setter修饰设置值的方法sett
成为群英会员,开启智能安全云计算之旅
立即注册关注或联系群英网络
7x24小时售前:400-678-4567
7x24小时售后:0668-2555666
24小时QQ客服
群英微信公众号
CNNIC域名投诉举报处理平台
服务电话:010-58813000
服务邮箱:service@cnnic.cn
投诉与建议:0668-2555555
Copyright © QY Network Company Ltd. All Rights Reserved. 2003-2020 群英 版权所有
增值电信经营许可证 : B1.B2-20140078 粤ICP备09006778号 域名注册商资质 粤 D3.1-20240008