xarray基础之计算篇
import xarray as xr
# 数据
f = '/home/mw/input/cru3546/cru_ts4.07.2021.2022.pre.dat.nc'
# 打开数据集
ds = xr.open_dataset(f)
dsWarning: ecCodes 2.21.0 or higher is recommended. You are running version 2.14.1<xarray.Dataset>
Dimensions: (lon: 720, lat: 360, time: 24)
Coordinates:
• lon (lon) float32 -179.8 -179.2 -178.8 -178.2 ... 178.8 179.2 179.8
• lat (lat) float32 -89.75 -89.25 -88.75 -88.25 ... 88.75 89.25 89.75
• time (time) datetime64[ns] 2021-01-16 2021-02-15 ... 2022-12-16
Data variables:
pre (time, lat, lon) float32 ...
stn (time, lat, lon) float64 ...
Attributes:
Conventions: CF-1.4
title: CRU TS4.07 Precipitation
institution: Data held at British Atmospheric Data Centre, RAL, UK.
source: Run ID = 2304141047. Data generated from:pre.2304141039.dtb
history: Fri 14 Apr 12:49:03 BST 2023 : User f098 : Program makegrid...
references: Information on the data is available at http://badc.nerc.ac...
comment: Access to these data is available to any registered CEDA user.
contact: support@ceda.ac.uk
平均
all_mean_pre = ds.pre.mean()
print('两年全球平均降水:', all_mean_pre)两年全球平均降水: <xarray.DataArray 'pre' ()>
array(61.810722, dtype=float32)first_mean_pre = ds.pre.sel(time=slice('2021-01-16', '2021-12-16')).mean()
print('2021年全球平均降水:', first_mean_pre)2021年全球平均降水: <xarray.DataArray 'pre' ()>
array(61.74274, dtype=float32)point_mean_pre = ds.pre.sel(time=slice('2021-01-16', '2021-12-16'),lat=32.75,lon=60.25).mean()
print('2021年某格点平均降水:', point_mean_pre)2021年某格点平均降水: <xarray.DataArray 'pre' ()>
array(12.158334, dtype=float32)
Coordinates:
lon float32 60.25
lat float32 32.75总和
sum_pre = ds.pre.sum()
print('两年全球总降水:', sum_pre)两年全球总降水: <xarray.DataArray 'pre' ()>
array(1.00014696e+08, dtype=float32)标准差
std_pre = ds.pre.std()
print('pre std:', std_pre)pre std: <xarray.DataArray 'pre' ()>
array(85.1408844)最大值与最大值坐标
max=ds.pre.max()
max<xarray.DataArray 'pre' ()>
array(2311.80004883)
argmax=ds.pre.argmax(dim=['lon','lat'])
argmax{'lon': <xarray.DataArray 'pre' (time: 24)>
array([612, 641, 641, 256, 204, 548, 509, 545, 205, 576, 520, 612, 716,
612, 612, 612, 556, 543, 506, 556, 641, 603, 205, 565])
Coordinates:
* time (time) datetime64[ns] 2021-01-16 2021-02-15 ... 2022-12-16,
'lat': <xarray.DataArray 'pre' (time: 24)>
array([199, 170, 170, 189, 185, 216, 207, 220, 188, 212, 205, 199, 144,
199, 199, 199, 204, 230, 218, 208, 168, 209, 188, 191])
Coordinates:
* time (time) datetime64[ns] 2021-01-16 2021-02-15 ... 2022-12-16}输出了二十四时间点的最大值坐标，下面选一个时间点
ds.pre[0].max()<xarray.DataArray 'pre' ()>
array(1128.09997559)
Coordinates:
time datetime64[ns] 2021-01-16
ds.pre[0].argmax(dim=['lat','lon']){'lat': <xarray.DataArray 'pre' ()>
array(199)
Coordinates:
time datetime64[ns] 2021-01-16,
'lon': <xarray.DataArray 'pre' ()>
array(612)
Coordinates:
time datetime64[ns] 2021-01-16}ds.pre[0][199,612]<xarray.DataArray 'pre' ()>
array(1128.1, dtype=float32)
Coordinates:
lon float32 126.2
lat float32 9.75
time datetime64[ns] 2021-01-16
Attributes:
long_name: precipitation
units: mm/month
correlation_decay_distance: 450.0
上面已找到最大值的经纬度，顺便画个图
import matplotlib.pyplot as plt
import numpy as np
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
from matplotlib.colors import LinearSegmentedColormap
import cmaps
lons=ds.lon
lats=ds.lat
data=ds.pre[0]
fig = plt.figure(figsize=(20,18))
# 设置地图投影为PlateCarree
ax = plt.axes(projection=ccrs.PlateCarree())
# 绘制填充等值线图
cs = ax.contourf(lons, lats, data, transform=ccrs.PlateCarree(), cmap=cmaps.radar, levels=np.linspace(data.min(), data.max(), 20))
# 添加网格线
ax.gridlines(draw_labels=True, dms=True, x_inline=False, y_inline=False)
# 添加海岸线
ax.coastlines()import matplotlib.pyplot as plt
import numpy as np
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
from matplotlib.colors import LinearSegmentedColormap
import cmaps
lons=ds.lon
lats=ds.lat
data=ds.pre[0]
fig = plt.figure(figsize=(20,10))
# 设置地图投影为PlateCarree
ax = plt.axes(projection=ccrs.PlateCarree())
# 绘制填充等值线图
cs = ax.contourf(lons, lats, data, transform=ccrs.PlateCarree(), cmap=cmaps.radar_1, levels=np.linspace(data.min(), data.max(), 20))
# 添加网格线
ax.gridlines(draw_labels=True, dms=True, x_inline=False, y_inline=False)
# 设置经纬度范围
ax.set_extent([120, 130, 5, 15], crs=ccrs.PlateCarree())
# 最大值打点
ax.scatter([data.lon[612]],[data.lat[199]],color='b',s=10)
# 添加海岸线
ax.coastlines()
ax.scatter([data.lon[612]],[data.lat[199]],color='b',s=20)
cbar = plt.colorbar(cs, orientation='vertical')/opt/conda/lib/python3.7/site-packages/cmaps/cmaps.py:3533: UserWarning: Trying to register the cmap 'radar_1' which already exists.
matplotlib.cm.register_cmap(name=cname, cmap=cmap)在东经126，北纬9的附近的小黑点就是最大值所在，格点数据放大看确实有点奇形怪状
还有对应的argmin，用法差不多就不多介绍了
获得下月减上月的数组
diff = ds.pre.diff(dim='time')
diff<xarray.DataArray 'pre' (time: 23, lat: 360, lon: 720)>
array([[[nan, nan, nan, ..., nan, nan, nan],
[nan, nan, nan, ..., nan, nan, nan],
[nan, nan, nan, ..., nan, nan, nan],
...,
[nan, nan, nan, ..., nan, nan, nan],
[nan, nan, nan, ..., nan, nan, nan],
[nan, nan, nan, ..., nan, nan, nan]],
[[nan, nan, nan, ..., nan, nan, nan],
[nan, nan, nan, ..., nan, nan, nan],
[nan, nan, nan, ..., nan, nan, nan],
...,
[nan, nan, nan, ..., nan, nan, nan],
[nan, nan, nan, ..., nan, nan, nan],
[nan, nan, nan, ..., nan, nan, nan]],
[[nan, nan, nan, ..., nan, nan, nan],
[nan, nan, nan, ..., nan, nan, nan],
[nan, nan, nan, ..., nan, nan, nan],
...,type=float32)Coordinates:
• lon (lon) float32 -179.8 -179.2 -178.8 -178.2 ... 178.8 179.2 179.8
• lat (lat) float32 -89.75 -89.25 -88.75 -88.25 ... 88.75 89.25 89.75
• time (time) datetime64[ns] 2021-02-15 2021-03-16 ... 2022-12-16
