如何用rioxarray绘制GeoTIFF并显示经纬度坐标？

问题背景 #

使用样本GeoTIFF文件绘图时，希望显示-180至180范围的经纬度坐标，而非当前的投影坐标系数值。现有代码及xarray.DataArray对象da的信息如下：

现有代码 #

import rioxarray
import matplotlib.pyplot as plt

fl = 'https://download.osgeo.org/geotiff/samples/gdal_eg/cea.tif'
da = rioxarray.open_rasterio(fl, masked=True)
da.plot()
plt.show()

da对象内容 #

<xarray.DataArray (band: 1, y: 515, x: 514)>
[264710 values with dtype=float32]
Coordinates:
  * band         (band) int64 1
  * x            (x) float64 -2.846e+04 -2.84e+04 ... 2.268e+03 2.328e+03
  * y            (y) float64 4.256e+06 4.256e+06 ... 4.225e+06 4.225e+06
    spatial_ref  int64 0
Attributes:
    AREA_OR_POINT:  Area
    scale_factor:   1.0
    add_offset:     0.0

解决方法 #

当前da的x、y是CEA投影坐标系下的数值，要显示经纬度，需将数据重投影到WGS84经纬度坐标系（EPSG:4326），具体步骤如下：

1. 重投影坐标系
   使用rioxarray的rio.reproject()方法将数据转换到EPSG:4326，转换后会自动生成对应经纬度的坐标值。

2. 绘图展示
   重投影后的数组绘图时，坐标轴会自动显示-180至180的经度、-90至90的纬度。

完整代码 #

import rioxarray
import matplotlib.pyplot as plt

fl = 'https://download.osgeo.org/geotiff/samples/gdal_eg/cea.tif'
da = rioxarray.open_rasterio(fl, masked=True)

# 重投影到WGS84经纬度坐标系
da_wgs84 = da.rio.reproject("EPSG:4326")

# 可选：将x/y重命名为lon/lat，让坐标含义更明确
da_wgs84 = da_wgs84.rename(x="lon", y="lat")

# 绘图，此时坐标轴显示经纬度
da_wgs84.plot()
plt.show()

说明 #

• 重投影过程中，rioxarray会基于原数据的空间参考信息（spatial_ref）自动完成坐标转换，无需手动计算。
• 转换后的da_wgs84的坐标范围会符合经纬度的标准区间，直接绘图即可显示目标坐标。

内容的提问来源于stack exchange，提问作者Rafael