我开始使用网格数据开始我的旅程,而我一直在尝试从带有cartopy的netcdf文件中绘制一些温度数据。我遵循了一些示例,但我不明白为什么我的地块中间有一条白线。 (我已经检查了数据,矩阵充满了数字,没有NaN)
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import numpy as np
import xarray as xr
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatterimport glob
data = xr.open_dataset('aux1.nc')
lat = data.lat
lon = data.lon
time = data.time
Temp = data.air
#Calculo la temperatura media anual
Tanual = Temp.resample(time="y").mean()
#Promedio de todos los meses
Tprom = Temp.mean(dim="time").values
#Grafico
fig = plt.figure(figsize=(10, 4))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())
ax.coastlines()
ax.set_global()
ct = ax.contourf(lon,lat,Tprom,transform=ccrs.PlateCarree(),cmap="bwr")
ax.gridlines()
cb = plt.colorbar(ct,orientation="vertical",extendrect='True')
cb.set_label("Temperatura [°C]")
ax.set_xticks(np.arange(-180,181,60), crs=ccrs.PlateCarree())
ax.set_yticks(np.arange(-90,91,30), crs=ccrs.PlateCarree())
lon_formatter = LongitudeFormatter(zero_direction_label=True)
lat_formatter = LatitudeFormatter()
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)
好问题!问题在于,对于大多数栅格化的气候数据而言,经度坐标类似于:
array([1.25, 3.75, 6.25, ..., 351.25, 353.75, 356.25, 358.75])
因此没有明确的longitude=0
点,这通常在绘图中为您提供一条细白线。我有时也会在已发表的论文(甚至是《自然》)中看到这个问题!
有很多方法可以解决此问题,但是最简单的方法是使用cartopy
程序包,该程序包具有一个名为add_cyclic_point
的实用程序,该实用程序基本上可以在longitude=0
点的任一侧插入数据。 (参考:https://scitools.org.uk/cartopy/docs/v0.15/cartopy/util/util.html)
此方法的唯一缺点是,当使用xarray
时,它意味着您必须手动提取数据,然后丢失元数据,因此,我编写了一个函数,以在保持元数据的同时保持其良好易用性。
from cartopy.util import add_cyclic_coord
import xarray as xr
def xr_add_cyclic_point(da):
"""
Inputs
da: xr.DataArray with dimensions (time,lat,lon)
"""
# Use add_cyclic_point to interpolate input data
lon_idx = da.dims.index('lon')
wrap_data, wrap_lon = add_cyclic_point(da.values, coord=da.lon, axis=lon_idx)
# Generate output DataArray with new data but same structure as input
outp_da = xr.DataArray(data=wrap_data,
coords = {'time': da.time, 'lat': da.lat, 'lon': wrap_lon},
dims=da.dims,
attrs=da.attrs)
return outp_da
示例
例如,如果我的初始DataArray看起来像:
<xarray.DataArray 'tas' (time: 60, lat: 90, lon: 144)>
[777600 values with dtype=float32]
Coordinates:
* lat (lat) float64 -89.49 -87.98 -85.96 -83.93 ... 85.96 87.98 89.49
* lon (lon) float64 1.25 3.75 6.25 8.75 11.25 ... 351.3 353.8 356.2 358.8
* time (time) object 1901-01-16 12:00:00 ... 1905-12-16 12:00:00
Attributes:
long_name: Near-Surface Air Temperature
units: K
valid_range: [100. 400.]
cell_methods: time: mean
standard_name: air_temperature
original_units: deg_k
original_name: t_ref
cell_measures: area: areacella
associated_files: baseURL: http://cmip-pcmdi.llnl.gov/CMIP5/dataLocation...
当我绘制时间平均值时,给出以下内容:
tas.mean(dim='time').plot.contourf()
现在,我可以使用我的函数来生成一个新的,内插的DataArray,如下所示:
wrapped_tas = xr_add_cyclic_point(tas)
wrapped_tas
<xarray.DataArray (time: 60, lat: 90, lon: 145)>
array([[[251.19466, 251.19469, 251.19472, ..., 251.19226, 251.19073,
251.19466], ...
[250.39403, 250.39468, 250.39961, ..., 250.39429, 250.39409,
250.39403]]], dtype=float32)
Coordinates:
* time (time) object 1901-01-16 12:00:00 ... 1905-12-16 12:00:00
* lat (lat) float64 -89.49 -87.98 -85.96 -83.93 ... 85.96 87.98 89.49
* lon (lon) float64 1.25 3.75 6.25 8.75 11.25 ... 353.8 356.2 358.8 361.2
Attributes:
long_name: Near-Surface Air Temperature
units: K
valid_range: [100. 400.]
cell_methods: time: mean
standard_name: air_temperature
original_units: deg_k
original_name: t_ref
cell_measures: area: areacella
associated_files: baseURL: http://cmip-pcmdi.llnl.gov/CMIP5/dataLocation...
您可以看到,经度坐标已从144-> 145的长度延长了一个点,这意味着它现在“环绕”了longitude=0
点。
此新的DataArray在绘制时给出了没有白线的图:)
wrapped_tas.mean(dim='time').plot.contour()
希望有帮助!! :)