关于GPM IMERG Final Half Hourly 数据问题

Question

我从官网上下载中国区的GPM半小时降水数据， nc格式。现使用python代码从nc数据中提取precipitation变量为tif数据，后面发现两个问题
1.用ArcMap查看方向是正确的，但用envi打开就会发生偏转。
2.我研究用python代码标出GPM中的降水位置，输入的是上面生成的正确的GPM.tif，但生成的可视化图还是翻转的原图降水区域明明集中在是左下角和右上角对角线
图一是原栅格图 很明显降水集中在左下角到右上角的对角线
图二是可视化栅格降水信息 降水集中到了左上角到右下角的对角线

下面是我用的代码：

import rasterio
import numpy as np
import matplotlib.pyplot as plt

# 1. 加载GPM数据
gpm_file = r"E:\XC\GPMjz_tif\3B-HHR.MS.MRG.3IMERG.20220515-S110000-E112959.0660.V07B.HDF5_clipped_precipitation_data.tif"
with rasterio.open(gpm_file) as gpm_src:
    gpm_data = gpm_src.read(1)  # 读取GPM数据
    gpm_meta = gpm_src.meta  # 获取GPM的元数据
    print("GPM Meta:", gpm_meta)  # 打印元数据，检查仿射变换矩阵

# 显示仿射变换矩阵
print("Affine Transform:", gpm_meta['transform'])

# 仿射变换矩阵一般是一个 3x2 矩阵：
# | a | b | x0 |
# | c | d | y0 |
# 其中：
# a, c = X轴方向的分辨率和旋转角度
# b, d = Y轴方向的分辨率和旋转角度
# x0, y0 = 参考点的原始坐标

# 2. 可视化原始GPM强降水区域
plt.imshow(gpm_data, cmap='Blues')
plt.title('Original GPM Strong Precipitation Mask')
plt.colorbar()
plt.show()

输出：

E:\AnacondaEnvs\XCpython\python.exe E:\project\pythonProject\对齐\经纬度顺序.py 
GPM Meta: {'driver': 'GTiff', 'dtype': 'float32', 'nodata': -9999.900390625, 'width': 52, 'height': 42, 'count': 1, 'crs': CRS.from_wkt('GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]'), 'transform': Affine(0.09999997008080576, 0.0, 117.99999543732288,
       0.0, 0.10000000930413967, 26.999999232408477)}
Affine Transform: | 0.10, 0.00, 118.00|
| 0.00, 0.10, 27.00|
| 0.00, 0.00, 1.00|
 
进程已结束，退出代码为 0
 

Answer 1

上午好☀️☀️☀️️
本答案参考通义千问

从你的描述来看，问题的核心在于数据的地理参考信息（Geo-referencing）与实际数据内容不一致，导致数据在不同的软件中显示时出现旋转或偏移现象。以下是问题分析及解决方案：

---

问题分析

1. ArcMap 和 ENVI 的行为差异

   • ArcMap 能正确解析 GeoTIFF 文件的元数据（包括仿射变换矩阵），因此能够正确显示数据的方向。
   • ENVI 在处理 GeoTIFF 文件时可能没有正确解析仿射变换矩阵，导致数据被错误地旋转或偏移。

2. Python 可视化结果翻转

   • 使用 rasterio 加载的数据中，仿射变换矩阵 (gpm_meta['transform']) 表明数据的 X 轴和 Y 轴方向的分辨率分别为正数（0.1°），并且原点位于左上角（x0=117.99999543732288, y0=26.999999232408477）。
   • 但是，可视化时默认的图像方向是从左上角开始，而数据的实际内容可能是从左下角开始。这导致了图像的翻转。

3. 数据存储顺序问题

   • rasterio.read() 默认返回的数据是按行优先存储的，而 GeoTIFF 文件的元数据可能表明数据存储顺序需要调整。

---

解决方案

1. 检查和修正仿射变换矩阵

仿射变换矩阵决定了数据的空间参考，确保其正确性是关键。可以通过以下方式验证和修正：

• 检查仿射变换矩阵的方向
  如果仿射变换矩阵的行列顺序不符合预期，可以手动调整矩阵。

• 修正矩阵
  如果矩阵的方向不正确，可以通过调整矩阵的参数来纠正。例如，将矩阵的行和列顺序互换。

from rasterio.transform import from_origin

# 原始仿射变换矩阵
affine_matrix = gpm_meta['transform']

# 创建一个新的仿射变换矩阵，确保数据方向正确
new_affine = from_origin(x_min=affine_matrix.c, y_max=affine_matrix.f, x_res=affine_matrix.a, y_res=-affine_matrix.e)

# 更新元数据
gpm_meta['transform'] = new_affine

---

2. 调整数据存储顺序

如果数据的实际存储顺序与可视化期望的顺序不同，可以通过 numpy.flip 或 numpy.transpose 来调整数据的排列。

# 调整数据的存储顺序
gpm_data = np.flip(gpm_data, axis=0)  # 将数据从左下角到右上角翻转为左上角到右下角

---

3. 重新保存 GeoTIFF 文件

在完成数据方向调整后，重新保存为 GeoTIFF 文件，确保后续软件能够正确解析。

# 保存修正后的 GeoTIFF 文件
with rasterio.open(
    "corrected_gpm.tif",
    "w",
    driver="GTiff",
    height=gpm_data.shape[0],
    width=gpm_data.shape[1],
    count=1,
    dtype=gpm_data.dtype,
    crs=gpm_meta['crs'],
    transform=new_affine,
    nodata=gpm_meta['nodata']
) as dst:
    dst.write(gpm_data, 1)

---

4. 可视化修正后的数据

最后，使用 matplotlib 或其他工具对修正后的数据进行可视化，确保方向正确。

# 可视化修正后的数据
plt.imshow(gpm_data, cmap='Blues')
plt.title('Corrected GPM Strong Precipitation Mask')
plt.colorbar()
plt.show()

---

完整代码示例

以下是完整的代码，包含所有步骤：

import rasterio
from rasterio.transform import from_origin
import numpy as np
import matplotlib.pyplot as plt

# 1. 加载GPM数据
gpm_file = r"E:\XC\GPMjz_tif\3B-HHR.MS.MRG.3IMERG.20220515-S110000-E112959.0660.V07B.HDF5_clipped_precipitation_data.tif"
with rasterio.open(gpm_file) as gpm_src:
    gpm_data = gpm_src.read(1)  # 读取GPM数据
    gpm_meta = gpm_src.meta  # 获取GPM的元数据
    print("GPM Meta:", gpm_meta)  # 打印元数据，检查仿射变换矩阵

# 检查仿射变换矩阵
print("Affine Transform:", gpm_meta['transform'])

# 修正仿射变换矩阵
new_affine = from_origin(
    x_min=gpm_meta['transform'].c,  # 左下角的经度
    y_max=gpm_meta['transform'].f,  # 左下角的纬度
    x_res=gpm_meta['transform'].a,  # 分辨率（经度方向）
    y_res=-gpm_meta['transform'].e  # 分辨率（纬度方向）
)
gpm_meta['transform'] = new_affine

# 调整数据存储顺序
gpm_data = np.flip(gpm_data, axis=0)  # 将数据从左下角到右上角翻转为左上角到右下角

# 重新保存 GeoTIFF 文件
with rasterio.open(
    "corrected_gpm.tif",
    "w",
    driver="GTiff",
    height=gpm_data.shape[0],
    width=gpm_data.shape[1],
    count=1,
    dtype=gpm_data.dtype,
    crs=gpm_meta['crs'],
    transform=new_affine,
    nodata=gpm_meta['nodata']
) as dst:
    dst.write(gpm_data, 1)

# 可视化修正后的数据
plt.imshow(gpm_data, cmap='Blues')
plt.title('Corrected GPM Strong Precipitation Mask')
plt.colorbar()
plt.show()

---

总结

通过上述步骤，可以解决 GPM 数据在不同软件中显示方向不一致的问题。关键点在于：

1. 检查并修正仿射变换矩阵，确保其正确反映数据的空间参考。
2. 调整数据存储顺序，使数据的排列符合可视化要求。
3. 重新保存 GeoTIFF 文件，确保后续软件能够正确解析。

如果仍有问题，请提供更多细节（如原始数据的具体范围、分辨率等），以便进一步分析！