Python:如何解决MODIS GRID(正弦投影/GCTP_SNSOID)的重投影问题?

2023-12-15 13:52:13

博客已同步微信公众号:GIS茄子;若博客出现纰漏或有更多问题交流欢迎关注GIS茄子,或者邮箱联系(推荐-见主页).

微信公众号

01 投影参数了解

这里我将用手上存有的MCD12Q1(土地利用数据集)和MOD11A2(地表温度数据集)说明隐藏在HDF4文件中的关键投影参数。

对于MODIS产品中的GRID数据集,一般都是使用USGS General Cartographic Transformation Package (GCTP,美国地质调查局通用制图转换包)的标准进行坐标系写入。

那么就需要说明GCTP对于各个坐标系的一些标准,例如LambertMercatorSinusoidal等,此处以Sinusoidal为例说明投影参数。

在说明之前,了解到MDOSI数据集(多HDF4文件)的投影参数一般可以通过软件查看(PanoplyHDFExplorer),当然也可以通过编程例如IDL的HDF4相关函数、Python中的xarraypyhdfgdal等都可以打开HDF4文件。

但是需要注意的是:Panoply在显示上会基于元数据(StructMetadata.0)进行一些数据集和属性的虚构,另外绘图也是如此(和实际的栅格矩阵或者数据集常规绘制出来有时不同甚至差异很大(例如自动纠正南北极显示<如源数据南极在上方>等)-这是因为Panoply在绘图上做了很多优化,这些优化多依靠地理信息)

如下:

Panoply软件正虚构了一些属性如上:Projection属性

而HDFExplorer软件则是正常显示:

HDFExplorer软件原汁原味的显示

所以通过编程手段获取时,我们实际上无法获取得到Projection这类虚构的属性,因为他们实际上不存在。(但是这不意味着HDFExplorer软件表现更好,对于NC文件该软件实际并不支持无法正常打开,功能也相对较少,现NASA似乎已经不进行该软件的更新和维护了)

好,我们回到问题,如何查看得到投影信息呢?既然是从元数据中进行虚构,那么就从元数据中寻找去,如下查看全局属性(StructMetadata.0):

StructMetadat属性

Projection表示当前数据集的投影坐标系是基于正弦投影(GCTP_SNSOID)建立的。

ProjParams表示投影参数:(6371007.181000,0,0,0,0,0,0,0,0,0,0,0,0);

但是上述投影参数是什么意思呢?

通过GCTP_SNSOID知道,这是GCTP标准下的正弦投影参数,因此需要查看关于GCTP的相关文档,这里给出相关链接:

这是关于GTCP下各个坐标系的投影参数详解:https://www.cmascenter.org/ioapi/documentation/all_versions/html/GCTP.html
PDF文件见:https://www.cmascenter.org/ioapi/documentation/all_versions/html/GCTP.pdf
NASA网站上也有提及:https://wiki.earthdata.nasa.gov/download/attachments/239406229/GeneralCartographicTransformationPackage_v2.pdf?version=1&modificationDate=1667483241998&api=v2
在网站中我们可以看到:

正弦投影参数说明述
可以发现:对于前面的ProjParams属性:(6371007.181000,0,0,0,0,0,0,0,0,0,0,0,0)知道,其中:
6371007.181000表示地球半径(下标为1处,下标从1开始);
投影中心的纬度为0°(下标为6,即赤道)
假东偏移量为0(下标为7, 即水平方向上添加的一个偏移值,用于调整坐标系统的原点)
假北偏移量为0(下标为8,即垂直方向上添加的一个偏移值,与假东偏移量类似)
其余下标的数值没有被使用或者说没有赋予实际含义,忽略即可。

但是,在MOD11A2数据集(地表温度数据集)中存在:

MOD11A2的元数据

在下标为9的位置为86400,目前尚不清楚原因,翻看第二版本的GCTP说明文档如下:

GeneralCartographicTransformationPackage_v2.pdf
虽然多了下标为5的中央经线,但是依旧没有关于第9位的参数,因此这里忽略该数值(实际重投影后与其它数据集进行对比,似乎不存在偏移和不确定现象)。
(如果对于对于第9位参数有想法或者知晓相关信息,请邮箱联系我或者微信公众号私信留言,这对我帮助很大,谢谢)

02 如何重投影?

基于给定的投影参数我建议通过proj4语法进行坐标系系统的创建,会比绝大多数方法更好且便于理解和操作,仅仅涉及简单的字符串操作。

例如前例子中的坐标系转换为proj4字符串为:

proj4_str = '+proj=sinu +lon_0=0.0 +lat_0=0.0 +x_0=0.0 +y_0=0.0 +R=6371007.181'

上述+proj表示投影类型,lon_0表示中央经线,lat_0表示中央纬线,x_0表示东偏,y_0表示北偏,R表示地球半径,更多信息查看Proj4官方文档。

时间关系,这里直接贴出代码:

def img_warp(src_img: np.ndarray, out_path: str, transform: list, src_proj4: str, out_res: float,
             dst_nodata: Union[int, float] = None, resample: str = 'nearest') -> None:
    """
    该函数用于对正弦投影下的栅格矩阵进行重投影(GLT校正), 得到WGS84坐标系下的栅格矩阵并输出为TIFF文件
    :param src_img: 待重投影的栅格矩阵
    :param out_path: 输出路径
    :param transform: 仿射变换参数([x_min, x_res, 0, y_max, 0, -y_res], 旋转参数为0是常规选项)
    :param out_res: 输出的分辨率(栅格方形)
    :param dst_nodata: 设置为NoData的数值
    :param out_type: 输出的数据类型
    :param resample: 重采样方法(默认是最近邻, ['nearest', 'bilinear', 'cubic'])
    :param src_proj4: 表达源数据集(src_img)的坐标系参数(以proj4字符串形式)
    :return: None
    """

    # 输出数据类型
    if np.issubdtype(src_img.dtype, np.integer):
        out_type = gdal.GDT_Int32
    elif np.issubdtype(src_img.dtype, np.floating):
        out_type = gdal.GDT_Float32
    else:
        raise ValueError("当前待校正数组类型为不支持的数据类型")
    resamples = {'nearest': gdal.GRA_NearestNeighbour, 'bilinear': gdal.GRA_Bilinear, 'cubic': gdal.GRA_Cubic}
    # 原始数据集创建(正弦投影)
    driver = gdal.GetDriverByName('MEM')  # 在内存中临时创建
    src_ds = driver.Create("", src_img.shape[1], src_img.shape[0], 1, out_type)  # 注意: 先传列数再传行数, 1表示单波段
    srs = osr.SpatialReference()
    srs.ImportFromProj4(src_proj4)  # 依据
    """
    对于src_proj4, 依据元数据StructMetadata.0知:
        Projection=GCTP_SNSOID; ProjParams=(6371007.181000,0,0,0,0,0,0,0,0,0,0,0,0)
    或数据集属性(MODIS_Grid_8Day_1km_LST/Data_Fields/Projection)知:
        :grid_mapping_name = "sinusoidal";
        :longitude_of_central_meridian = 0.0; // double
        :earth_radius = 6371007.181; // double
    """
    src_ds.SetProjection(srs.ExportToWkt())  # 设置投影信息
    src_ds.SetGeoTransform(transform)  # 设置仿射参数
    src_ds.GetRasterBand(1).WriteArray(src_img)  # 写入数据
    # 重投影信息(WGS84)
    dst_srs = osr.SpatialReference()
    dst_srs.ImportFromEPSG(4326)
    # 重投影
    dst_ds = gdal.Warp(out_path, src_ds, dstSRS=dst_srs, xRes=out_res, yRes=out_res, dstNodata=dst_nodata,
                       outputType=out_type, multithread=True, format='GTiff', resampleAlg=resamples[resample])
    if dst_ds:  # 释放缓存和资源
        dst_ds.FlushCache()
        src_ds, dst_ds = None, None

上述过程还不包括如何提取元数据的信息,这里仅贴出代码,至于更详细的内容正在更新中。

from pyhdf.SD import SD
hdf = SD(hdf_path)
# 获取元数据
metadata = hdf.__getattr__('StructMetadata.0')
# 获取角点信息
ul_pt= [float(x) for x in re.findall(r'UpperLeftPointMtrs=\((.*)\)', metadata)[0].split(',')]
lr_pt = [float(x) for x in re.findall(r'LowerRightMtrs=\((.*)\)', metadata)[0].split(',')]
x_mins.append(ul_pt[0])
x_maxs.append(lr_pt[0])
y_mins.append(lr_pt[1])
y_maxs.append(ul_pt[1])
# 计算分辨率
col = int(re.findall(r'XDim=(.*?)\n', metadata)[0])
row = int(re.findall(r'YDim=(.*?)\n', metadata)[0])
x_res = (lr_pt[0] - ul_pt[0]) / col
y_res = (ul_pt[1] - lr_pt[1]) / row
# 获取数据集的坐标系参数并转化为proj4字符串格式
projection_param = [float(_param) for _param in re.findall(r'ProjParams=\((.*?)\)', metadata)[0].split(',')]
mosaic_img_proj4 = "+proj={} +R={:0.4f} +lon_0={:0.4f} +lat_0={:0.4f} +x_0={:0.4f} " \
                   "+y_0={:0.4f} ".format('sinu', projection_param[0], projection_param[4], projection_param[5],
                                          projection_param[6], projection_param[7])
# 关闭文件, 释放资源
hdf.end()

简单修改了一点,可能存在纰漏,如果留有问题请与我联系。

文章来源:https://blog.csdn.net/m0_63001937/article/details/135012787
本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。