准备工作
Hello 写这篇博文主要是最近在数据处理上需要快速的分段统计栅格数据中的数值分布，所以就此分享给大家。
为了方便大家实地测试，我使用的数据都是开源获取的方便大家自己测试哈！！！
数据来源
气温数据
气温数据是老演员了，这次还是使用的NCEP/NCAR的数据气温数据作为本次教程的测试数据，不过下载NCEP气温数据是NC格式的，我这里使用ArcGis转换了一下导出为tif格式的。NC文件下载地址：点我
；如果大家不是很会使用ArcGis的转化的话，我这里提供了一个处理好的供大家测试：点我。这样设置应该是免费的，如果不是免费的话大家戳我一下我给大家放网盘下载。
代码
其实想要对栅格数据进行分组统计并不是很难的。最核心的思想就是把tif转为数组，然后再把数组降维成一行就行了。
需要用的库
# 加载GDAL
from osgeo import gdal
import numpy as np
import pandas as pd
# 画图用的库
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
import matplotlib.pyplot as plt
import cartopy.feature as cfeature
from cartopy.util import add_cyclic_point
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
获取从TIF中获取经纬度
tif_path =r'E:\CSDN\air_layer_ra1.tif' # 文件目录
dataset = gdal.Open(tif_path) # 打开tif
adfGeoTransform = dataset.GetGeoTransform() # 读取地理信息
print(adfGeoTransform[0]) # 左上角x坐标(经度)
print(adfGeoTransform[3]) # 左上角y坐标(纬度)
nXSize = dataset.RasterXSize # 列数
nYSize = dataset.RasterYSize # 行数
print(nXSize)
print(nYSize)
# 生成经纬度列表
# 警告下面这种方法无法生成带有投影坐标系的经纬度列表
# 因为投影坐标系涉及到椭球体
# 我这个是WGS84坐标系的所以可以直接这么生成
# 生成经纬度
lat = []
lon = []
for i in range(nXSize):
lon.append(adfGeoTransform[0] + (i+1) * adfGeoTransform[1])
for i in range(nYSize):
lat.append(adfGeoTransform[3] + i * adfGeoTransform[5])
lon = np.array(lon)
lat = np.array(lat)
数据过滤
air_temp = dataset.ReadAsArray() # 读取tif中的数据
air_temp = air_temp.astype(np.float32) # 设定数据类型
air_temp[air_temp < -9999] = np.nan # 将无效值替换为nan ArcGis中默认的无效值是一个非常小的负数，过滤掉
air_temp = np.array(air_temp)
air_max = np.nanmax(air_temp) # 获取最大值
air_min = np.nanmin(air_temp) # 获取最小值
print(air_temp.shape)
print(air_max,air_min)
TIF数据可视化
extent=[-180,180,-90,90]#限定绘图范围
projection = ccrs.PlateCarree()
# 设置画布大小
fig = plt.figure(figsize=(10,8),dpi=300)
ax = plt.axes(projection=projection)
air_cyc, lon_cyc = add_cyclic_point(air_temp, coord = lon)
ax.set_extent(extent, crs=ccrs.PlateCarree()) # 设置绘图范围
ax.xaxis.set_major_formatter(LongitudeFormatter(zero_direction_label=True))
ax.yaxis.set_major_formatter(LatitudeFormatter())
#画线范围和间隔
ax.set_xticks(np.arange(-180,180.1,60),crs=projection) # 设置x轴刻度
ax.set_yticks(np.arange(-90,90.1,30),crs=projection) # 设置y轴刻度
c = ax.contourf(lon_cyc,lat,air_cyc,cmap='coolwarm',transform=ccrs.PlateCarree()) # 绘制等值线
ax.add_feature(cfeature.COASTLINE.with_scale('50m'), linewidth=0.6)
plt.colorbar(c,ax=ax,orientation='horizontal',shrink=0.8,pad=0.05,extend='both')
plt.title('Air_Temp',fontsize=16)
数值分布
air_temp_list = air_temp.reshape(1,-1)[0].tolist() #转化成列表
air_temp_list = [x for x in air_temp_list if np.isnan(x) == False] # 过滤掉NAN
air_temp_list_num = len(air_temp_list) #获取数据个数
# 设定横坐标轴的最大值和最小值范围
if int(air_max/10)*10==air_max:
air_max_10 = air_max
else:
air_max_10 = int(air_max/10)*10+10
air_min_10 = int(air_min/10)*10
air_level = np.arange(air_min_10,air_max_10+0.1,10)
print(len(air_temp_list))
air_temp_level_sta = pd.cut(air_temp_list, air_level) # 核心步骤，统计每个分组的数据个数
print(air_temp_level_sta.value_counts())
分组可视化
x = []
y = []
for i in range(len(air_temp_level_sta.value_counts())):
x.append(str(air_temp_level_sta.value_counts().index[i]))
y.append(air_temp_level_sta.value_counts().values[i])
plt.figure(figsize=(8,6),dpi=300)
plt.plot(x,y,label='分组统计',marker='o')
plt.xticks(rotation=45)
# 下面这一组代码是为了在图上显示统计信息
# 用不上的话可以直接注释掉
plt.text(0.5,0.5,'最大值：'+str(air_max),transform=plt.gca().transAxes)
plt.text(0.5,0.4,'最小值：'+str(air_min),transform=plt.gca().transAxes)
plt.text(0.5,0.3,'平均值：'+str(np.mean(air_temp_list)),transform=plt.gca().transAxes)
plt.text(0.5,0.2,'标准差：'+str(np.std(air_temp_list)),transform=plt.gca().transAxes)
plt.text(0.5,0.1,'总数：'+str(air_temp_list_num),transform=plt.gca().transAxes)
plt.legend()
完整代码奉上
# 加载GDAL
from osgeo import gdal
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
import matplotlib.pyplot as plt
import cartopy.feature as cfeature
from cartopy.util import add_cyclic_point
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
tif_path =r'E:\CSDN\air_layer_ra1.tif'
dataset = gdal.Open(tif_path) # 打开tif
adfGeoTransform = dataset.GetGeoTransform() # 读取地理信息
print(adfGeoTransform[0]) # 左上角x坐标(经度)
print(adfGeoTransform[3]) # 左上角y坐标(纬度)
nXSize = dataset.RasterXSize # 列数
nYSize = dataset.RasterYSize # 行数
print(nXSize)
print(nYSize)
# 生成经纬度列表
# 警告下面这种方法无法生成带有投影坐标系的经纬度列表
# 因为投影坐标系涉及到椭球体
# 我这个是WGS84坐标系的所以可以直接这么生成
lat = []
lon = []
for i in range(nXSize):
lon.append(adfGeoTransform[0] + (i+1) * adfGeoTransform[1])
for i in range(nYSize):
lat.append(adfGeoTransform[3] + i * adfGeoTransform[5])
lon = np.array(lon)
lat = np.array(lat)
air_temp = dataset.ReadAsArray()
air_temp = air_temp.astype(np.float32)
air_temp[air_temp < -9999] = np.nan # 将无效值替换为nan ArcGis中默认的无效值是一个非常小的负数，过滤掉
air_temp = np.array(air_temp)
air_max = np.nanmax(air_temp)
air_min = np.nanmin(air_temp)
print(air_temp.shape)
print(air_max,air_min)
extent=[-180,180,-90,90]#限定绘图范围
projection = ccrs.PlateCarree()
# 设置画布大小
fig = plt.figure(figsize=(10,8),dpi=300)
ax = plt.axes(projection=projection)
air_cyc, lon_cyc = add_cyclic_point(air_temp, coord = lon)
ax.set_extent(extent, crs=ccrs.PlateCarree()) # 设置绘图范围
ax.xaxis.set_major_formatter(LongitudeFormatter(zero_direction_label=True))
ax.yaxis.set_major_formatter(LatitudeFormatter())
#画线范围和间隔
ax.set_xticks(np.arange(-180,180.1,60),crs=projection) # 设置x轴刻度
ax.set_yticks(np.arange(-90,90.1,30),crs=projection) # 设置y轴刻度
c = ax.contourf(lon_cyc,lat,air_cyc,cmap='coolwarm',transform=ccrs.PlateCarree()) # 绘制等值线
ax.add_feature(cfeature.COASTLINE.with_scale('50m'), linewidth=0.6)
plt.colorbar(c,ax=ax,orientation='horizontal',shrink=0.8,pad=0.05,extend='both')
plt.title('Air_Temp',fontsize=16)
air_temp_list = air_temp.reshape(1,-1)[0].tolist()
air_temp_list = [x for x in air_temp_list if np.isnan(x) == False]
air_temp_list_num = len(air_temp_list)
if int(air_max/10)*10==air_max:
air_max_10 = air_max
else:
air_max_10 = int(air_max/10)*10+10
air_min_10 = int(air_min/10)*10
air_level = np.arange(air_min_10,air_max_10+0.1,10)
print(len(air_temp_list))
air_temp_level_sta = pd.cut(air_temp_list, air_level)
print(air_temp_level_sta.value_counts())
x = []
y = []
for i in range(len(air_temp_level_sta.value_counts())):
x.append(str(air_temp_level_sta.value_counts().index[i]))
y.append(air_temp_level_sta.value_counts().values[i])
plt.figure(figsize=(8,6),dpi=300)
plt.plot(x,y,label='分组统计',marker='o')
plt.xticks(rotation=45)
# 下面这一组代码是为了在图上显示统计信息
# 用不上的话可以直接注释掉
plt.text(0.5,0.5,'最大值：'+str(air_max),transform=plt.gca().transAxes)
plt.text(0.5,0.4,'最小值：'+str(air_min),transform=plt.gca().transAxes)
plt.text(0.5,0.3,'平均值：'+str(np.mean(air_temp_list)),transform=plt.gca().transAxes)
plt.text(0.5,0.2,'标准差：'+str(np.std(air_temp_list)),transform=plt.gca().transAxes)
plt.text(0.5,0.1,'总数：'+str(air_temp_list_num),transform=plt.gca().transAxes)
plt.legend()
扫一扫
1953
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
