Python-气象-大气科学-可视化绘图系列(一)——利用xarray读取netCDF文件并画图(代码+示例)

本文原创链接:https:www.cnblogs.com/zhanling/p/12192978.html

import numpy as np
import xarray as xr
import cartopy.crs as ccrs
import cartopy.feature as cfeat
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import matplotlib.pyplot as plt


ds = xr.open_dataset('2039071310.003.nc')
t = ds['value']
lons = ds.lon.data
lats = ds.lat.data
temp = xr.DataArray(t.data.T, coords=[lats,lons], dims=['latitude','longitude'])

# 创建画图空间
proj = ccrs.PlateCarree() #创建投影
fig = plt.figure(figsize=(16,9)) #创建页面
ax = fig.subplots(1, 1, subplot_kw={'projection': proj}) #子图
# 设置地图属性:加载国界、海岸线、河流、湖泊
ax.add_feature(cfeat.BORDERS.with_scale('50m'), linewidth=0.8, zorder=1)
ax.add_feature(cfeat.COASTLINE.with_scale('50m'), linewidth=0.6, zorder=1)
ax.add_feature(cfeat.RIVERS.with_scale('50m'), zorder=1)
ax.add_feature(cfeat.LAKES.with_scale('50m'), zorder=1)
# 设置网格点属性
gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
linewidth=1.2, color='k', alpha=0.5, linestyle='--')
gl.xlabels_top = False #关闭顶端标签
gl.ylabels_right = False #关闭右侧标签
gl.xformatter = LONGITUDE_FORMATTER #x轴设为经度格式
gl.yformatter = LATITUDE_FORMATTER #y轴设为纬度格式
# 设置colorbar
cbar_kwargs = {
'orientation': 'horizontal',
'label': 'Potential',
'shrink': 0.8,
}
# 画图
levels = np.arange(0,1,0.1)
temp.plot.contourf(ax=ax, levels=levels, cmap='Spectral_r',
cbar_kwargs=cbar_kwargs, transform=ccrs.PlateCarree())
plt.savefig('test.jpg')


---------------------
 

示例效果(强对流概率预报结果)
---------------------
著作权归作者所有。
商业转载请联系作者获得授权,非商业转载请注明出处。
作者:斩羚_HermionX

 

python详细编写:有从北纬54-56经度81到83,空间长度为0.1度,从1994年到2023年7月1到31号0-23时的24个温度数据。数据保存在txt文档。坐标点大小是21*21,(54,81)对应的坐标是(1,1),(54,81.1)对应的坐标是(1,2),(54,81.2)对应的坐标是(1,3),...,(54,83)对应的坐标是(1,21),...,(54.1,81)对应的是(2,1),(54.1,81.1)对应的是(2,2),(54.1,81.2)对应的是(2,3),...,(54.1,83)对应的是(2,21),以此类推(56,81)对应的是(21,1),(56,81.1)对应的是(21,2),(56,81.2)对应的是(21,3),...,(56,83)对应的是(21,21)。第步建向量组ξ:在坐标点(1,1)的数据是(54,81)对应的温度,ξ1=(ξ0,1,ξ1,1,ξ2,1,ξ3,1...ξ23,1),ξ2=(ξ0,2,ξ1,2,ξ2,2,ξ3,2...ξ23,2),ξ3=(ξ0,3,ξ1,3,ξ2,3,ξ3,3...ξ23,3)...ξ31=(ξ0,31,ξ1,31,ξ2,31,ξ3,31...ξ23,31),ξ=(ξ1,ξ2,ξ3...ξ31)这里的ξ大小为31*21, ξ1,ξ2...ξ31是每年7月份1- 31号,里面包含的0-23的24个元素是每隔小时的气温。以此类推,这样的坐标点有21*21个,k,l=1,2,..21。再代入1994-2023年的数据。第二步求加权方差的计算:M(ξ-)*(ξ-)^T=∑,这里∑的大小为31*31是托普利兹矩阵。第三步求均值 (1/30)*∑n=1,...,30*((1/31)*∑i=1,...,31*ξi,n(k,l))=Mξ。k,l=1,..,21,这里k,l表示坐标点。(1/30)*∑n=1,...,30*((1/31)*∑i=1,...,31*ξi,n(k,l)) 表示双重求和的平均值。从坐标(1,1)直到坐标(21,21)的ξ,对从1到31天进行求和,除以31来计算平均值。对n从1到30(1994-2023共30年)求和除以30来计算平均值。Mξ=(μ...μ)总共31个μ。μ=(μ1,μ2...μ24).Mξ大小是744。坐标点(k,l)(k,l)为Mξ在k,l坐标点的值。第四步M(ξ-)*(ξ-)^T= Mξ(ξ^T) - Mξ*M(ξ^T )。求Mξ(ξ^T)((ξ^T))(h)=(1/30)*∑n=1,...,30*((1/(31-h))*∑i=1,...,(31-h)*ξi,n*(ξi+h,n)^T),h=0,1,...,30。(ξi,n)这里的i表示天数,1到31。n表示年数,从1到30(1994-2023共30年).ξi+h,n这里的i+h也表示天数,n表示年数 ,h为0到30。Mξ(ξ^T) 表示矩阵Mξ与矩阵ξ的转置的乘积.第五步求(1/30)*∑n=1,...,30*((1/31)*∑i=1,...,31*ξi,n验证(1/30)*∑n=1...30*((1/31)*∑i=1,...,31*ξi,n是否等于μ,μ=(μ1,μ2...μ24).第六步求R.R=(R0 R1 R2...R30; R1^T, R0 R1...R29;...;R30^T R29^T...R0),其中R0=((ξ^T))(h=0)-μ*(μ^T),R1=((ξ^T))(h=1)-μ*(μ^T),R2=((ξ^T))(h=2)-μ*(μ^T)...Rk=((ξ^T))(h=k)-μ*(μ^T).R是个31x31的矩阵,由R0到R30(大小为24x24)组成.Mξ(ξ^T)(h) 表示矩阵Mξ与矩阵ξ的转置的乘积,乘积中的元素是由数据进行求和和平均值计算得到。第六步重新给ξ添加上下标的值,上标是i,n,下标是k,I.i是1~31号,n是1~30年.k,l是坐标点(21*21)(kl)=(1/30)*∑n=1...30*((1/31)*∑i=1...31*ξi,n,k,l). A=(1/30)*∑n=1...30*((1/31)*∑i=1...31*ξi,n,k1,l1*ξi,n,k2,l2)-μk1,l1*μk2,l2.其中ξk1,l1,ξk2,l2是每天中的任意两个时间点的温度。每天任选两个时间点的温度,总共31(31天)*30()对。第七步画图。在21*21的坐标中中点是(11,11),找到由中点到坐标中任意点的相关系数画图呈现出来。
最新发布
03-11
<think>嗯,用户想用Python处理1994到2023年北纬54-56度、东经81-83度区域的逐小时温度数据。首先,我需要确定数据处理的步骤。根据问题,需要创建温度向量组ξ,计算加权方差,求日平均温度,协方差矩阵,以及绘制相关系数图。 首先,数据读取和预处理是关键。用户的数据可能是CSV或NetCDF格式,比如气象常用的NetCDF。需要用到xarray库来高效处理这些多维数据。然后筛选指定经纬度范围和时间范围,这里可以用xarray的sel方法结合切片操作。 接下来,处理缺失值。气象数据常有缺失,比如传感器故障或传输问题。需要用fillna或插值方法处理,但用户的问题里没有特别提到,可能需要先检查数据完整性。 然后,构建温度数据向量组ξ。这里可能需要将每个坐标点的温度按时间序列排列,形成个向量。可能需要reshape操作,把数据转换成(时间点,空间点)的二维数组。 计算加权方差∑。用户提到加权方差,但需要明确权重是什么。通常可能是时间或空间上的权重。假设是时间上的,比如不同时间的权重不同,可以用numpy的加权函数,例如np.average配合weights参数,然后计算方差。 日平均温度μ需要将逐小时数据按天聚合。可以用xarray的resample方法,按天重采样,取平均值。这会得到每个坐标点每天的平均温度。 协方差矩阵R的计算,需要各个坐标点温度数据之间的协方差。使用numpy的cov函数,但要注意数据的排列方式,通常是每行个变量,每列个观测值,所以可能需要转置输入矩阵。 绘图部分,基于中心点(11,11)和其他点的相关性。需要提取中心点的温度序列,计算与其他各点的相关系数,可以用np.corrcoef。然后用matplotlib绘制热图或散点图,显示相关性。可能需要调整坐标轴标签,显示实际的经纬度值。 另外,用户提到的引用中有提到数据质量的问题,比如GPS偏移。虽然这里是温度数据,但也需要注意数据的准确性,比如异常值的处理。可能需要在预处理步骤中加入数据清洗,比如去除超出合理范围的温度值。 还有,用户引用的案例中使用了串口接收GPS数据,但这里处理的是历史数据,可能不需要实时采集,但数据存储和读取方式需要注意。比如,处理大量逐小时数据时,要考虑内存管理,可能分块处理。 最后,代码的整合和测试。确保每步的输出符合预期,尤其是维度转换和重采样步骤。可能需要在Jupyter Notebook中逐步运行,检查中间结果是否正确。</think>根据需求,分步骤实现如下(假设数据已存储为NetCDF格式): ### 、数据预处理 1. **数据加载与筛选** ```python import xarray as xr # 读取数据集 ds = xr.open_dataset('temperature_data.nc') # 筛选区域:北纬54-56,东经81-83 data = ds.sel(latitude=slice(54,56), longitude=slice(81,83)) # 提取温度变量(假设变量名为't2m') temp_data = data['t2m'] ``` 2. **构建温度向量组$ξ$** ```python # 将三维数据(时间,纬度,经度)转为二维矩阵(时间,空间点) ξ = temp_data.stack(space=['latitude','longitude']).transpose('time','space') ``` ### 二、统计分析 3. **计算加权方差$Σ$** ```python import numpy as np # 定义时间权重(示例:指数衰减权重) weights = np.exp(-0.1 * np.arange(len(ξ.time))) weights /= weights.sum() # 计算加权方差 Σ = np.average((ξ - ξ.mean(dim='time'))**2, axis=0, weights=weights) ``` 4. **计算日平均温度$μ$** ```python # 按日重采样求平均 daily_mean = ξ.resample(time='1D').mean() # 提取具体坐标点(假设坐标为(11,11)) μ = daily_mean.sel(space=(54.5,81.5)) # 根据实际坐标索引调整 ``` ### 三、协方差矩阵计算 5. **构建协方差矩阵$R$** ```python # 转置为(空间点,时间) R = np.cov(ξ.values.T) ``` ### 四、可视化实现 6. **相关系数图绘制** ```python import matplotlib.pyplot as plt # 计算相关系数矩阵 corr_matrix = np.corrcoef(ξ.values.T) # 提取中心点(11,11)的相关系数 center_idx = 11 * ξ.space.shape[1] + 11 # 假设网格排列方式 correlations = corr_matrix[center_idx] # 绘制热力图 plt.figure(figsize=(12,8)) plt.imshow(correlations.reshape(ξ.latitude.size, ξ.longitude.size), cmap='coolwarm', extent=[81,83,54,56]) plt.colorbar(label='Correlation Coefficient') plt.xlabel('Longitude') plt.ylabel('Latitude') plt.title('Temperature Correlation with Center Point (54.5N,81.5E)') plt.show() ``` ### 关键问题说明: 1. **数据预处理**需注意时间连续性,如遇数据中断应进行插值处理[^1] 2. **权重设计**需根据具体需求选择,气象领域常用余弦纬度加权 3. **坐标索引转换**需通过`xarray`的网格坐标系统实现精确定位 4. **可视化精度**可通过`cartopy`或`folium`提升地理信息展示效果
评论 20
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值