处理GEE输出的tif图像

需求分析：S5P的官方数据，极难处理。借助GEE平台，下载其中的数据（tif格式）并处理为通用的nc文件

数据：在gee中export的tif文件。其在gee中的可视化：

---

（2022/11/23更新）

  喵了个咪的，发现xarray直接可以打开.tif文件，并进行可视化绘图，以2019-01-10的GEE MCD19A2为例：

# （1）使用xarray包下的open_rasterio读取（没有则conda install rasterio）

Dir2 = 'F:\\2_Joint_Estimation\\Source Data\\GEE_AOD2019'
files2 = os.listdir(Dir2)
df = xr.open_rasterio(os.path.join(Dir2, files2[0]))
t = df[0]*0.001
df2 = xr.open_rasterio(os.path.join(Dir2, files2[1]))
t2 = df2[0]*0.001
t_com = xr.concat([t, t2], dim='x') # 按经度拼接

t_com.plot(cmap="Spectral_r", vmin=0, vmax=1)  # Spectral的互逆配色
# plt.savefig('gee-aod_' + t_com.time.data.astype(str)[0][:10] + '.png')
plt.show()

# (2) 或者使用更新包 rioxarray (conda install rioxarray -c conda-forge)
"""使用rioxarray打开文件无warning，但没有concat拼接功能"""
# df = rioxarray.open_rasterio(os.path.join(Dir2, files2[0]))

（nc VS tif）

  左侧为使用原答案方法将tif转为nc、绘图；右侧为直接使用tif读取绘图。即左侧空缺值使用nan填充，右侧使用0填充。

一、将tif转为nc文件

from osgeo import gdal
import xarray as xr
import matplotlib.pyplot as plt


def tiff2nc(path):
    data = gdal.Open("temp_file/CHINA_NO2_1Km_201906.tif")
    im_width = data.RasterXSize  # 获取宽度，数组第二维，左右方向元素长度，代表经度范围
    im_height = data.RasterYSize  # 获取高度，数组第一维，上下方向元素长度，代表纬度范围
    im_bands = data.RasterCount  # 波段数
    """
    GeoTransform 的含义：
        影像左上角横坐标：im_geotrans[0]，对应经度
        影像左上角纵坐标：im_geotrans[3]，对应纬度
    
        遥感图像的水平空间分辨率(纬度间隔)：im_geotrans[5]
        遥感图像的垂直空间分辨率(经度间隔)：im_geotrans[1]
        通常水平和垂直分辨率相等
    
        如果遥感影像方向没有发生旋转，即上北下南，则 im_geotrans[2] 与 im_geotrans[4] 为 0
    
    计算图像地理坐标：
        若图像中某一点的行数和列数分别为 row 和 column，则该点的地理坐标为：
            经度：xGeo = im_geotrans[0] + col * im_geotrans[1] + row * im_geotrans[2]
            纬度：yGeo = im_geotrans[3] + col * im_geotrans[4] + row * im_geotrans[5]
    """
    im_geotrans = data.GetGeoTransform()  # 获取仿射矩阵，含有 6 个元素的元组

    im_proj = data.GetProjection()  # 获取地理信息

    """
    GetRasterBand(bandNum)，选择要读取的波段数，bandNum 从 1 开始
    ReadAsArray(xoff, yoff, xsize, ysize)，一般就按照下面这么写，偏移量都是 0 ，返回 ndarray 数组
    """
    im_data = data.GetRasterBand(1).ReadAsArray(xoff=0, yoff=0, win_xsize=im_width, win_ysize=im_height)
    # 根据im_proj得到图像的经纬度信息
    im_lon = [im_geotrans[0] + i * im_geotrans[1] for i in range(im_width)]
    im_lat = [im_geotrans[3] + i * im_geotrans[5] for i in range(im_height)]

    im_nc = xr.DataArray(im_data, coords=[im_lat, im_lon], dims=['lat', 'lon'])
    return im_nc

day_nc = tiff2nc('temp_file/CHINA_NO2_1Km_201906.tif')

# 使用xarray自带的绘图函数进行可视化
day_nc.plot()
plt.show()

二、转化后nc文件的可视化：

三、总结：
  两幅图的效果不能说一模一样，也能说毫无差别 [允悲]。这里估摸着可能是colorbar的原因，gee调色板：palette: [‘black’, ‘blue’, ‘purple’, ‘cyan’, ‘green’, ‘yellow’, ‘red’]。
  存在的不足：（1）gee导出图像的时间略长，该样例导出时间约8min/张；（2）导出的图像暂未找到时间信息，需要手动添加。。
参考文章：『GDAL』读写TIFF文件