【问题标题】:plotting data from netcdf with cartopy isnt plotting data at 0 longitude使用 cartopy 绘制来自 netcdf 的数据不是在 0 经度处绘制数据
【发布时间】:2020-07-09 12:10:05
【问题描述】:

我从网格数据开始我的旅程,我一直在尝试从带有 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)

https://i.stack.imgur.com/aPHz4.jpg

【问题讨论】:

    标签: python netcdf python-xarray cartopy


    【解决方案1】:

    好问题!问题在于,对于大多数网格化气候数据,经度坐标类似于:

    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()
    

    希望有帮助!! :)

    【讨论】:

      猜你喜欢
      • 1970-01-01
      • 2021-05-28
      • 2019-02-07
      • 2013-12-29
      • 2014-08-13
      • 2020-06-11
      • 1970-01-01
      • 1970-01-01
      • 2016-06-09
      相关资源
      最近更新 更多