Python:如何解决MODIS GRID(正弦投影/GCTP_SNSOID)的重投影问题?
博客已同步微信公众号:GIS茄子
;若博客出现纰漏或有更多问题交流欢迎关注GIS茄子
,或者邮箱联系(推荐-见主页).
01 投影参数了解
这里我将用手上存有的MCD12Q1
(土地利用数据集)和MOD11A2
(地表温度数据集)说明隐藏在HDF4文件中的关键投影参数。
对于MODIS产品中的GRID数据集,一般都是使用USGS General Cartographic Transformation Package (GCTP
,美国地质调查局通用制图转换包)的标准进行坐标系写入。
那么就需要说明GCTP对于各个坐标系的一些标准,例如Lambert
,Mercator
,Sinusoidal
等,此处以Sinusoidal
为例说明投影参数。
在说明之前,了解到MDOSI数据集(多HDF4
文件)的投影参数一般可以通过软件查看(Panoply
、HDFExplorer
),当然也可以通过编程例如IDL的HDF4相关函数、Python中的xarray
、pyhdf
、gdal
等都可以打开HDF4文件。
但是需要注意的是:Panoply在显示上会基于元数据(StructMetadata.0
)进行一些数据集和属性的虚构,另外绘图也是如此(和实际的栅格矩阵或者数据集常规绘制出来有时不同甚至差异很大(例如自动纠正南北极显示<如源数据南极在上方>等)-这是因为Panoply在绘图上做了很多优化,这些优化多依靠地理信息)
如下:
而HDFExplorer软件则是正常显示:
所以通过编程手段获取时,我们实际上无法获取得到Projection这类虚构的属性,因为他们实际上不存在。(但是这不意味着HDFExplorer软件表现更好,对于NC文件该软件实际并不支持无法正常打开,功能也相对较少,现NASA似乎已经不进行该软件的更新和维护了)
好,我们回到问题,如何查看得到投影信息呢?既然是从元数据中进行虚构,那么就从元数据中寻找去,如下查看全局属性(StructMetadata.0
):
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数据集(地表温度数据集)中存在:
在下标为9的位置为86400,目前尚不清楚原因,翻看第二版本的GCTP说明文档如下:
虽然多了下标为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()
简单修改了一点,可能存在纰漏,如果留有问题请与我联系。
本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。 如若内容造成侵权/违法违规/事实不符,请联系我的编程经验分享网邮箱:veading@qq.com进行投诉反馈,一经查实,立即删除!