本文主要是介绍ENVI IDL：如何MODIS GRID产品进行批量镶嵌、重投影(GLT校正)？，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！
00 前言
最近比赛繁忙，催更MODIS相关处理的渐多，恰逢最近关于MODIS GRID产品的相关处理还没开始做，今天做个了断。
在处理之前我们先了解MODIS GRID产品的数据格式和使用方法等等相关内容。
0.1 数据格式
我们首先理解一下文件名各个部分（以第一个标红文件示例）的含义：
这是文件名：MOD11B3.A2021335.h31v07.061.2022004015800.hdf
MOD11B3是数据集的名称，一般通过该数据集你可以很轻松的知道关于它的一些信息（通过google、bing等），这是NASA对于它的信息说明：https://ladsweb.modaps.eosdis.nasa.gov/missions-and-measurements/products/MOD11B3/。那么就仅该部分而言我们可以获得的信息有哪些？MOD表示安装在Terra卫星上的MODIS(Moderate Resolution Imaging Spectroradiometer，中等分辨率成像光谱仪)传感器，类似的MYD表示安装在Aqua卫星上的MODIS传感器，而MCD表示二者的综合；至于11指代的Land Surface Temperature地表温度（显然这是一个Level-3级产品）；B3似乎表示数据集为月度产品，在NASA官网的介绍中：Each LST&E pixel value in the MOD11B3 is a simple average of all the corresponding values from the MOD11B1 collected during the month period. 整体理解，MOD11B3表示MODIS Terra卫星的月度陆地表面温度和发射率合成产品（包括白天和黑夜）。
A20121335表示数据集上的原始数据拍摄的时间，即2012年第335天。时间格式为YYYYDDD。
h31v07表示数据集在GRID格网上的图块编码位置（此处又特质SIN正弦投影GRID格网上），此即第31列第7行的图块。关于全球各个图块的相对位置见下图：，关于更详细的信息查看NASA官网：https://modis-land.gsfc.nasa.gov/MODLAND_grid.html，其中还包括Lambert GRID、Linear GRID等GRID格网的相关说明。
061表示算法版本，一般而言选择最新版本算法更好精度更高。此处为6.1版本。
2022004015700表示产品生产的时间，也就是你手上拿到的这儿HDF文件生成的时间，它不是卫星拍摄影像的时间。因为这是一个三级产品，由卫星拍摄影像到你这个三级产品有一系列的算法进行生产。此处表示2022年第4天凌晨1点57分00秒生产的，时间格式为YYYYDDDHHMMSS上述所有时间均遵循UTC世界协调时，确保所有产品都在一个全球统一的时间上进行，如果想要转化为中国北京时间应为UTC + 8时，因为数据集为月度数据，时间分辨率不算高，因此后续处理并没有考虑这一步骤。
.hdf这是文件后缀，一般.hdf表示HDF4文件，.h5或者hdf5等表示HDF5文件，但并非绝对。
如果想要更清楚了解MOD1B3数据集的数据格式和使用说明，请下载用户指南：
https://lpdaac.usgs.gov/documents/715/MOD11_User_Guide_V61.pdf
0.2 如何查看数据内容
关于数据内容，通过Panoply软件进行HDF文件的打开和显示。下载官网为：https://www.giss.nasa.gov/tools/panoply/download/。现已经更新到5.2.10版本。
但是需要注意其需要java11以上环境，需要自行去Oracle Java下载java11以上环境。
或者使用HDFView软件，详情查看：https://www.hdfgroup.org/downloads/hdfview/。
或者使用HDFExp软件，但是NASA似乎已经停用了，可以自行在网络上寻找。
0.3 数据内容说明
这是关于MOD11B3内部层次结构的详细说明：https://ladsweb.modaps.eosdis.nasa.gov/filespec/MODIS/6/MOD11B3
01 思路
1.1 思路1
- 对于每一个HDF4文件提取目标数据集，并对该数据集基于地理定位信息（通过
StructMetadata.0元数据集中查找）进行GLT校正(即重投影为WGS84坐标系)并输出为TIFF文件临时存储； - 接着对所有TIFF文件进行
镶嵌（即均值处理）；
但是实际处理后发现存在镶嵌缝：
这种缝隙产生的原因查看思路3；
思路1解决该方法就是采用最简单的滑动窗口填充，注意不是最近邻、双线性、三次卷积等插值，仅仅是一个3✖3的滑动窗口对所有像素进行遍历，对于缺失的像素使用以该像元为中心的3✖3的窗口内其他有效像元值的简单平均值作为该缺失像元的填充值。
1.2 思路2
深究原理？如下为中间数据：
由此可见，不是镶嵌导致的问题，而是在多个镶嵌图像镶嵌之前就已经不能恰好无缝衔接。这说明生产各个镶嵌图像之前就存在一些纰漏。于是我们尝试去查找是否为原始影像上就存在不衔接的情况：
两张影像一个为h25v04(在上)一个为h25v05(在下)，是上下衔接的两张影像，查看角点信息发现，h25v05的左上角点的Y坐标为4447802.079066与h25v04的右下角点的Y坐标完全一致。说明本身数据不存在问题，而是我们进行处理时存在一些纰漏。
那么在镶嵌之前我们进行的操作就是将HDF4文件中的目标数据集基于角点信息进行GLT校正并输出为WGS84坐标系的GeoTIFF文件，这说明这一步存在问题。具体为什么存在问题查看思路3；
思路2的解决办法就是我们先在SIN(正弦投影)坐标系下进行镶嵌，镶嵌成一张大的栅格矩阵，然后基于这个大栅格矩阵基于角点信息进行GLT校正。即将原先的先校正后镶嵌换为先镶嵌后校正。
1.3 思路3
对于原始数据正常接边，但是在完成GLT校正之后存在缝隙，实际上是没有认真考虑角点信息的含义。
这是一份角点信息的示例：
普遍认为左上角点和右下角点指代的为下方所示位置：
那我们姑且认为这种说法是正确的。
那么在实际求解每一像元的坐标时，又如此操作：
for ix=0, ds_size[0] - 1 do x_prj_coor[ix, *] = ul_coor[0] + x_res * ix
for ix=0, ds_size[1] - 1 do y_prj_coor[*, ix] = ul_coor[1] - y_res * ix
那么如此得到的每一像素点坐标实际上是该像元的左上角点位置的坐标，而非该像元的中心处坐标，若想要得到像元中心处坐标，应：
for ix=0, ds_size[0] - 1 do x_prj_coor[ix, *] = ul_coor[0] + x_res * ix + x_res / 2.0 ; 表示像元的中心位置
for ix=0, ds_size[1] - 1 do y_prj_coor[*, ix] = ul_coor[1] - y_res * ix - y_res / 2.0
那么当进行GLT校正时，poly_2d函数存在参数pixel_center，官方说明为：
PIXEL_CENTER
Set this keyword to a floating-point value giving the offset of the center of each input pixel. The default is 0.0 which indicates the lower-left corner of the pixel. A typical value would be 0.5, which would be the center of the pixel.
即该参数用于给定每一个输入像元的偏置量，默认偏置量是0.0，其表示像元的左下角点位置；一个比较经典的数值是0.5，它表示像元的中心位置。
按照上述说明，但是似乎一个浮点型数值pixel_center无法表示像元的左上角位置，它似乎是从数学坐标系的左下角开始，从X和Y轴均增加pixel_center(一个像素的百分之)，0.5恰好是中心位置。
如果我们进行GLT校正之前对于坐标计算是像元中心位置，那么GLT校正时pixel_center就应该设置为0.5，如下：
target_warped = poly_2d(target, k_col, k_row, interp_code[interp], $col_count_out, row_count_out, missing=!values.F_NAN, pixel_center=0.5)
不进行上述操作得到结果如下：
进行上述操作结果如下：
但是似乎还是存在一些拼接痕迹，那么回到前面我们交流的内容:普遍认为左上角点和右下角点指代的为下方所示位置那我们姑且认为这种说法是正确的。
我认为这种说法是错误的，这是一些参考资料：
这是HDF-EOS官网论坛下关于角点信息的激烈探讨：https://hdfeos.org/forums/showthread.php?t=460；
这是NASA下MODLand关于角点的简单说明：https://landweb.modaps.eosdis.nasa.gov/cgi-bin/developer/tilemap.cgi；
但是总的来说，应该可以这么理解，对于HDF文件的元数据中，应该存在两个参数（但往往存在一个因为非此即彼）：GridOrigin和PixelRegistration。
参数PixelRegistration有两个可选，一是HDFE_CENTER，二是HDFE_CORNER。如果指定HDFE_CENTER，那么说明角点信息表示左上角像元的中心坐标；如果指定HDFE_CORNER，那么还需要指定GridOrigin参数，它包括HDFE_GD_UL、HDFE_GD_UR、HDFE_GD_LL、HDFE_GD_LR依次为左上、右上、左下、右下位置；
但是一般指定了GridOrigin参数就可能不会指定PixelRegistration参数为HDFE_CORNER。
例如在OMI Aura数据中：
例如在MODIS11B3数据中：
因此上述角点信息实际上是表示：
但是，在代码中如何体现呢？
是的，仅仅修改这几行代码，另外计算坐标时注意为中心像元，最后GLT校正的pixel_center参数进行设置为0.5.整个过程为先GLT校正后镶嵌即可无缝隙，既不需要插值，也不需要调换顺序先镶嵌后GLT校正。
为什么这样可以呢？如果按照上述理解，那么相邻俩文件之间其实存在一个重叠行和重叠列，因此不会产生缝隙。
对于思路3基本上在思路1的基础上更换几行代码即可，因此无代码展示。自行理解。
03 代码
3.1 先GLT校正后镶嵌再填充
对于填充函数时间有限，为了更好避免对边界进行扩边，还可以采取阈值例如窗口内至少存在5个以上有效值才可以进行填充等等，一定程度上避免对边界进行扩边膨胀，另外窗口大小尽量小，3×3即可。太大也会导致边界膨胀。
另外这部分代码难度比较小，加之有现成代码，所以基本上无注释，见谅。另外前期写代码没有规范，所以一些地方尾巴处理不是很好，不够简洁。另外对一开始的时间提取采用的自定义类进行，就很pythonic，文末我会贴出自定义的date类。不喜欢这种面向对象的时间处理方法可以自行修改为函数加循环即可。
modis_grid_mosaic：主程序
; @Author : ChaoQiezi
; @Time : 2023年10月23日-下午9:43:49
; @Email : chaoqiezi.one@qq.com; 该程序用于 ···pro modis_grid_mosaic; 准备工作in_dir = 'D:\Objects\JuniorFallTerm\IDLProgram\Experiments\ExperimentalData\chapter_3\MOD11B3\'out_dir = 'D:\Objects\JuniorFallTerm\IDLProgram\Experiments\ExperimentalData\chapter_3\IDL_out_midterm_method1\'if ~file_test(out_dir) then file_mkdir, out_dirtarget_name = 'LST_Day_6km'sf_name = 'scale_factor'add_name = 'add_offset'range_name = 'valid_range'out_res = 0.06 ; 度; 检索img_paths = file_search(in_dir, '*.hdf', count=img_count)img_dates = strmid(file_basename(img_paths), 9, 7)for ix=0, img_count - 1 do img_dates[ix] = (date(ydays=img_dates[ix])).strftime('%Y-%M')img_dates_unique = img_dates.uniq()foreach img_date, img_dates_unique do beginprint, '正在镶嵌: ' + img_datemosaic_ixs = where(img_dates eq img_date); 循环-几何校正 out_warped_paths = out_dir + file_basename(img_paths[mosaic_ixs], '.hdf') + '.tiff' foreach mosaic_ix, mosaic_ixs, out_ix do begincur_path = img_paths[mosaic_ix]modis_grid_geo_transform, cur_path, cur_lon, cur_lat; 获取目标数据集target = read_h4(cur_path, target_name, /double)scale_factor = read_h4(cur_path, target_name, attr_name=sf_name)add_offset = read_h4(cur_path, target_name, attr_name=add_name)valid_range = read_h4(cur_path, target_name, attr_name=range_name); 数据集预处理invalid_ixs = where((target lt valid_range[0]) or (target gt valid_range[1]))target[invalid_ixs] = !values.F_NANtarget = target * scale_factor[0] + add_offset[0]; 几何校正.img_warp, target, cur_lon, cur_lat, out_res, target_warped, geo_info; 输出write_tiff, out_warped_paths[out_ix], target_warped, geotiff=geo_info, /floatprint, file_basename(out_warped_paths[out_ix], '.tiff'), '--几何校正并输出TIFF影像完毕'endforeach; 镶嵌img_mosaic, out_warped_paths, mosaic_img, geo_infofile_delete, out_warped_paths ; 删除中间文件; 插值window_interp, mosaic_img, mosaic_interp_img, window_size=3; 输出out_path = out_dir + img_date + '.tiff'write_tiff, out_path, mosaic_interp_img, geotiff=geo_info, /floatendforeach
end
modis_grid_geo_transform：读取MODIS GRID产品的角点信息并返回经纬度数据集
;+
; 函数用途:
; 该程序用于读取MODIS GRID产品的角点信息并返回经纬度数据集
; 函数参数:
; h4_path: MODIS GRID文件的路径
; extract_lon: 返回的经度数据集
; extract_lat: 返回的纬度数据集
; ds_size: (关键字参数)包含列号、行号的二元数组, 如果没有传入从文件中读取
;-
pro modis_grid_geo_transform, h4_path, extract_lon, extract_lat, ds_size=ds_size ; 获取元数据h4_id = hdf_sd_start(h4_path, /read)metadata_ix = hdf_sd_attrfind(h4_id, 'StructMetadata.0')hdf_sd_attrinfo, h4_id, metadata_ix, data=metadata; 获取行列号if ~keyword_set(ds_size) then beginds_size = make_array(2, /double)ix_one = strpos(metadata, 'XDim')ix_two = strpos(metadata, 'YDim')ix_three = strpos(metadata, 'UpperLeftPointMtrs')ds_size[0] = (strsplit(strmid(metadata, ix_one, ix_two - ix_one), '=', /extract))[1]ds_size[1] = (strsplit(strmid(metadata, ix_two, ix_three - ix_two), '=', /extract))[1]endif else beginds_size = double(ds_size) ; 转换为双精度, 避免后续超出范围endelse; 获取角点信息ix_one = strpos(metadata, 'UpperLeftPointMtrs')ix_two = strpos(metadata, 'LowerRightMtrs')ix_three = strpos(metadata, 'Projection')ul_coor = strsplit(strmid(metadata, ix_one, ix_two - ix_one), '=(,)', /extract)dr_coor = strsplit(strmid(metadata, ix_two, ix_three - ix_two), '=(,)', /extract)ul_coor = double(ul_coor[1:2])dr_coor = double(dr_coor[1:2])x_res = (dr_coor[0] - ul_coor[0]) / ds_size[0]y_res = (ul_coor[1] - dr_coor[1]) / ds_size[1]
; x_res = (dr_coor[0] - ul_coor[0]) / (ds_size[0] - 1) ; 对于思路3可更换此
; y_res = (ul_coor[1] - dr_coor[1]) / (ds_size[1] - 1)
; dr_coor[0] += x_res
; dr_coor[1] -= y_res; 计算-投影坐标数据集x_prj_coor = make_array(ds_size, /double)y_prj_coor = make_array(ds_size, /double)for ix=0, ds_size[0] - 1 do x_prj_coor[ix, *] = ul_coor[0] + x_res * ix + x_res / 2.0 ; 表示像元的中心位置for ix=0, ds_size[1] - 1 do y_prj_coor[*, ix] = ul_coor[1] - y_res * ix - y_res / 2.0
; for ix=0, ds_size[0] - 1 do x_prj_coor[ix, *] = ul_coor[0] + x_res * ix ; 老师写, 与poly_2d-pixel_center存在矛盾
; for ix=0, ds_size[1] - 1 do y_prj_coor[*, ix] = ul_coor[1] - y_res * ixx_prj_coor = reform(x_prj_coor, ds_size[0] * ds_size[1], /overwrite) ; /overwrite避免复制加快速度y_prj_coor = reform(y_prj_coor, ds_size[0] * ds_size[1], /overwrite); 进行reform是为了适配map_proj_inverse函数的输入: An n-element vector containing the x values. ; 投影坐标数据集 ==> 经纬度数据集prj_info = map_proj_init('Sinusoidal', /gctp, sphere_radius=6371007.181, $center_longitude=0.0, false_easting=0.0, false_northing=0.0); 需要说明的是gctp表示U.S. Geological Survey's General Cartographic Transformation Package (GCTP)来进行投影而不是IDL自己的投影库geo_coor = map_proj_inverse(x_prj_coor, y_prj_coor, map_structure=prj_info)extract_lon = geo_coor[0, *].reform(ds_size)extract_lat = geo_coor[1, *].reform(ds_size)
end
img_warp： 基于经纬度数据集对目标数据集进行GLT校正(重投影-WGS84)
;+
; 函数用途:
; 基于经纬度数据集对目标数据集进行几何校正(重投影-WGS84)
; 函数参数:
; target: 待校正的目标数据集
; lon: 对应目标数据集的经度数据集
; lat: 对应目标数据集的纬度数据集
; out_res: 输出的分辨率(°)
; target_warped: 校正后的目标数据集
; geo_info: 校正后的目标数据集对应的地理结构体
; degree<关键字参数: 5>: 多项式的次数
; interp<关键字参数: nearest>: 插值算法(包括: nearest, linear, cublic)
; sub_percent<关键字参数: 0.1>: 默认使用10%的点位进行几何校正
;-
pro img_warp, target, lon, lat, out_res, target_warped, geo_info, degree=degree, interp=interp, $sub_percent=sub_percent; 获取基本信息ds_size = size(target)lon_lat_corner = hash($ ; 考虑到min()max()为角点像元的中心处而非四个角点的边界'min_lon', min(lon) - out_res / 2.0d, $'max_lon', max(lon) + out_res / 2.0d, $'min_lat', min(lat) - out_res / 2.0d, $'max_lat', max(lat) + out_res / 2.0d)col_count_out = ceil((lon_lat_corner['max_lon'] - lon_lat_corner['min_lon']) / out_res)row_count_out = ceil((lon_lat_corner['max_lat'] - lon_lat_corner['min_lat']) / out_res)interp_code = hash($'nearest', 0, $'linear', 1, $'cublic', 2)if ~keyword_set(interp) then interp='nearest'if ~keyword_set(degree) then degree=5if ~keyword_set(sub_percent) then sub_percent = 0.1; 原始的行列号矩阵row_ori = make_array(ds_size[1:2], /integer)col_ori = make_array(ds_size[1:2], /integer)for ix=0, ds_size[1]-1 do col_ori[ix, *] = ixfor ix=0, ds_size[2]-1 do row_ori[*, ix] = ix; 校正后的行列号矩阵col_warp = floor((lon - lon_lat_corner['min_lon']) / out_res)row_warp = floor((lon_lat_corner['max_lat'] - lat) / out_res); 获取sub_percent的均匀采样点索引all_count = ds_size[1] * ds_size[2]sample_count = floor(all_count * sub_percent)sample_ix = floor((findgen(sample_count) / double(sample_count)) * all_count)polywarp, col_ori[sample_ix], row_ori[sample_ix], col_warp[sample_ix], row_warp[sample_ix], $degree, k_col, k_rowtarget_warped = poly_2d(target, k_col, k_row, interp_code[interp], $col_count_out, row_count_out, missing=!values.F_NAN, pixel_center=0.5)geo_info={$MODELPIXELSCALETAG: [out_res, out_res, 0.0], $ ; 分辨率MODELTIEPOINTTAG: [0.0, 0.0, 0.0, lon_lat_corner['min_lon'], lon_lat_corner['max_lat'], 0.0], $ ; 角点信息GTMODELTYPEGEOKEY: 2, $ ; 设置为地理坐标系GTRASTERTYPEGEOKEY: 1, $ ; 像素的表示类型, 北上图像(North-Up)GEOGRAPHICTYPEGEOKEY: 4326, $ ; 地理坐标系为WGS84GEOGCITATIONGEOKEY: 'GCS_WGS_1984', $GEOGANGULARUNITSGEOKEY: 9102} ; 单位为度
end
img_mosaic: 对多个文件基于地理位置进行镶嵌/均值运算
;+
; 函数用途:
; 该函数用于对多个文件基于地理位置进行镶嵌/均值运算
; 函数参数:
; img_paths: 所有需要镶嵌的影像的绝对路径(字符串数组形式)
; mosaic_img: 输出的镶嵌好的影像
; geo_info: 镶嵌好的影像的地理结构体(包含地理定位等信息)
;-
pro img_mosaic, img_paths, mosaic_img, geo_info; 准备lon_lat_set = {$min_lon: list(), $max_lat: list(), $min_lat: list(), $max_lon: list()}mosaic_imgs = list(); 循环获取镶嵌的统一空间范围foreach img_path, img_paths do beginimg = read_tiff(img_path, geotiff=geo_info)img_res = geo_info.MODELPIXELSCALETAGimg_tiepoint = geo_info.MODELTIEPOINTTAGimg_size = size(img)lon_lat_set.min_lon.add, img_tiepoint[3]lon_lat_set.max_lat.add, img_tiepoint[4]lon_lat_set.max_lon.add, (img_tiepoint[3]) + img_size[1] * img_res[0]lon_lat_set.min_lat.add, (img_tiepoint[4]) - img_size[2] * img_res[1]endforeachmin_lon = min(lon_lat_set.min_lon.toarray())max_lat = max(lon_lat_set.max_lat.toarray())min_lat = min(lon_lat_set.min_lat.toarray())max_lon = max(lon_lat_set.max_lon.toarray()); 镶嵌col_count = ceil((max_lon - min_lon) / img_res[0])row_count = ceil((max_lat - min_lat) / img_res[1])foreach img_path, img_paths do beginimg = read_tiff(img_path, geotiff=geo_info)img_tiepoint = geo_info.MODELTIEPOINTTAGimg_size = size(img)cur_mosaic_img = make_array(col_count, row_count, /double, value=!values.F_NAN) start_col = floor((img_tiepoint[3] - min_lon) / img_res[0])start_row = floor((max_lat - img_tiepoint[4]) / img_res[1])end_col = start_col + img_size[1] - 1end_row = start_row + img_size[2] - 1cur_mosaic_img[start_col: end_col, start_row: end_row] = imgmosaic_imgs.add, cur_mosaic_imgendforeachmosaic_imgs = mosaic_imgs.toarray(); 计算-输出mosaic_img = mean(mosaic_imgs, dimension=1, /nan) ; 可能警告Floating illegal operand, 原因为dimension=1的维度上求取均值遇到全为nan的情况if check_math(mask=128) eq 128 then beginprint, '警告: Floating illegal operand, 原因为镶嵌时dimension=1的维度上求取均值遇到全为nan的情况'endifgeo_info={$MODELPIXELSCALETAG: [img_res[0], img_res[1], 0.0], $ ; 分辨率MODELTIEPOINTTAG: [0.0, 0.0, 0.0, min_lon, max_lat, 0.0], $ ; 角点信息GTMODELTYPEGEOKEY: 2, $ ; 设置为地理坐标系GTRASTERTYPEGEOKEY: 1, $ ; 像素的表示类型, 北上图像(North-Up)GEOGRAPHICTYPEGEOKEY: 4326, $ ; 地理坐标系为WGS84GEOGCITATIONGEOKEY: 'GCS_WGS_1984', $GEOGANGULARUNITSGEOKEY: 9102} ; 单位为度
end
window_interp: 对栅格矩阵中缺失值基于滑动窗口进行填充
;+
; 函数用途:
; 该函数用于对栅格矩阵中缺失值基于滑动窗口进行填充
; 函数参数:
; dataset: 需要进行填补的栅格矩阵
; dataset_interp: 输出的经填补好的栅格矩阵
; window_size(默认: 3): 窗口大小(奇数)
;-
pro window_interp, dataset, dataset_interp, window_size=window_sizeds_size = size(dataset, /dimensions)ds_type = size(dataset, /type); 边界填充if ~keyword_set(window_size) then window_size = 3padding_size = window_size / 2ds_size += padding_size * 2dataset_pad = padding(dataset, padding_size)dataset_interp = padding(dataset, padding_size); 插值for col_ix=padding_size, ds_size[0] - padding_size - 1 do beginfor row_ix=padding_size, ds_size[1] - padding_size - 1 do begin; 若不是NAN跳过if ~finite(dataset_pad[col_ix, row_ix], /nan) then continue; 对NAN进行插值window_ds = dataset_pad[col_ix-padding_size: col_ix+padding_size, $row_ix-padding_size: row_ix+padding_size]valid_ix = where(~finite(window_ds, /nan), count)if (count eq 0) then continuedataset_interp[col_ix, row_ix] = mean(window_ds, /nan)endforendfor; no paddingdataset_interp = dataset_interp[padding_size:(ds_size[0] - padding_size - 1), $padding_size:(ds_size[1] - padding_size - 1)]
end
3.2 先镶嵌再GLT校正
这一部分是之前写，但是今天又整理了一下。需要注意的是，除了最开始对于时间处理使用到了自定义时间类外（需要自己修改，或者加上我后续的自定义时间类代码），其他地方均从底层开始写，不再封装任何函数，因此逻辑上看起来比较琐屑。
当然，代码中还有很多需要注意的点，因为顺序颠倒不仅仅是将代码上下简单换一下，还需要一些逻辑自洽，甚至注意的点比_3.1_部分相当，这里时间关系不再说明。
; @Author : ChaoQiezi
; @Time : 2023年11月6日-上午8:39:51
; @Email : chaoqiezi.one@qq.com; 该程序用于 ···(不套用任何自定义函数实现,仅可看帮助完成, 不可查看其他代码)pro modis_grid_mosaicstart_time = systime(1); 准备工作in_dir = 'D:\Objects\JuniorFallTerm\IDLProgram\Experiments\ExperimentalData\chapter_3\MOD11B3\' ; 输入路径out_dir = 'D:\Objects\JuniorFallTerm\IDLProgram\Experiments\ExperimentalData\chapter_3\IDL_out_midterm_method3\' ; 输出路径if ~file_test(out_dir) then file_mkdir, out_dir ; 如果out_dir输出文件夹不存在那么创建target_name = 'LST_Day_6km' ; 镶嵌的目标数据集 - 名称sf_name = 'scale_factor' ; 目标数据集的比例系数 scale factor - 名称add_name = 'add_offset' ; 目标数据集的偏置量add offset - 名称range_name = 'valid_range' ; 目标数据集的有效范围 - 名称out_res = 0.055d ; 输出分辨率(约等于6KM)degree = 5 ; 默认五次多项式(GLT校正)interp_code = 0 ; 0: 最近邻插值, 1: 线性插值, 2: 三次卷积插值sub_percent = 0.1 ; 使用控制点数量百分比(用于GLT校正)raw_rows = 200 ; 单-目标数据集的行列号raw_cols = 200img_paths = file_search(in_dir, '*.hdf') ; 检索所有HDF文件img_dates = list(strmid(file_basename(img_paths), 9, 7), /extract) ; 获取所有文件的日期并转化为列表形式存储img_dates= img_dates.map(lambda(a: (date(ydays=a)).strftime('%Y-%M'))) ; 利用自定义类date进行日期转化, 使用匿名函数+map对所有日期操作img_dates_unique = (img_dates.toarray()).uniq() ; 对日期唯一化; 对相同日期的影像依次进行镶嵌、GLT校正最终输出为GeoTIFF文件foreach img_date, img_dates_unique do beginmosaic_ixs = where(img_dates eq img_date, /null) ; 在当前循环日期下的所有下标mosaic_paths = img_paths[mosaic_ixs] ; 获取当前循环日期下的所有路径; 用于存储当前日期下所有镶嵌的图像的角点信息x_mins = list()x_maxs = list()y_mins = list()y_maxs = list(); 获取大镶嵌图像的角点信息和分辨率foreach mosaic_path, mosaic_paths do begincur_time = systime(1); 读取h4_id = hdf_sd_start(mosaic_path, /read)metadata_ix = hdf_sd_attrfind(h4_id, 'StructMetadata.0')hdf_sd_attrinfo, h4_id, metadata_ix, data=metadatahdf_sd_end, h4_id; 获取角点信息one_pos = strpos(metadata, 'UpperLeftPointMtrs=')two_pos = strpos(metadata, 'LowerRightMtrs=')three_pos = strpos(metadata, 'Projection=GCTP_SNSOID')ul_point = strmid(metadata, one_pos, two_pos - one_pos)lr_point = strmid(metadata, two_pos, three_pos - two_pos)ul_point = strsplit(ul_point, '=(,)', /extract)lr_point = strsplit(lr_point, '=(,)', /extract)x_mins.add, double(ul_point[1])x_maxs.add, double(lr_point[1])y_mins.add, double(lr_point[2])y_maxs.add, double(ul_point[2])endforeach
; raw_x_res = (x_maxs[0] - x_mins[0]) / raw_cols ; 老师所想
; raw_y_res = (y_maxs[0] - y_mins[0]) / raw_rows
; x_min = min(x_mins.toarray())
; x_max = max(x_maxs.toarray())
; y_min = min(y_mins.toarray())
; y_max = max(y_maxs.toarray())raw_x_res = (x_maxs[0] - x_mins[0]) / (raw_cols - 1) ; 严格的角点信息raw_y_res = (y_maxs[0] - y_mins[0]) / (raw_rows - 1)x_min = min(x_mins.toarray())x_max = max(x_maxs.toarray()) + raw_x_resy_min = min(y_mins.toarray()) - raw_y_resy_max = max(y_maxs.toarray())mosaic_cols = ceil((x_max - x_min) / raw_x_res)mosaic_rows = ceil((y_max - y_min) / raw_y_res); 对当前循环日期下的每一影像进行数据集和属性的读取 ==> 为镶嵌做准备mosaic_imgs = list()foreach mosaic_path, mosaic_paths, mosaic_ix do begin; 读取h4_id = hdf_sd_start(mosaic_path, /read)ds_ix = hdf_sd_nametoindex(h4_id, target_name)ds_id = hdf_sd_select(h4_id, ds_ix)hdf_sd_getdata, ds_id, targetsf_ix = hdf_sd_attrfind(ds_id, sf_name)add_ix = hdf_sd_attrfind(ds_id, add_name)range_ix = hdf_sd_attrfind(ds_id, range_name)hdf_sd_attrinfo, ds_id, sf_ix, data=scale_factorhdf_sd_attrinfo, ds_id, add_ix, data=add_offsethdf_sd_attrinfo, ds_id, range_ix, data=valid_rangehdf_sd_endaccess, ds_idhdf_sd_end, h4_id; 数据集预处理target = double(target) ; 避免后续赋值NAN出现浮点数警告invalid_ixs = where((target lt valid_range[0]) or (target gt valid_range[1])) ; 提取无效值下标target[invalid_ixs] = !values.F_NAN ; 赋值NANtarget = target * scale_factor[0] + add_offset[0] ; 换算; 将数据添加mosaic_img = make_array(mosaic_cols, mosaic_rows, /double, value=!values.F_NAN)x_min_cur = x_mins[mosaic_ix]x_max_cur = x_maxs[mosaic_ix]y_min_cur = y_mins[mosaic_ix]y_max_cur = y_maxs[mosaic_ix]x_start = floor((x_min_cur - x_min) / raw_x_res)y_start = floor((y_max - y_max_cur) / raw_y_res)x_end = x_start + raw_cols - 1y_end = y_start + raw_rows - 1mosaic_img[x_start:x_end, y_start:y_end] = targetmosaic_imgs.add, mosaic_imgendforeachmosaic_imgs = mean(mosaic_imgs.toarray(), dimension=1, /nan) ; 镶嵌(实为镶嵌) ==> 此处占用时间最长
; ; 监测是否存在浮点数非法警告(mask=128代号表示该警告) ==> % Program caused arithmetic error: Floating illegal operand
; if check_math(mask=128) eq 128 then begin
; print, '深呼吸, 某一像元上均为NAN导致求取该像元均值出现浮点数非法为正常现象'
; endif; 获取镶嵌图像的经纬度数据集(部分-有效); 为避免镶嵌图像上的NAN(一部分为确实是数值缺失, 另一部分可能是无该坐标在坐标系上, 理解参考FY4A宇宙部分像素), 此处仅有效像元部分进行GLT校正valid_pixel_ix = where(~finite(mosaic_imgs, /nan)) ; 获取有效像元的下标xs = make_array(mosaic_cols, mosaic_rows, /double)ys = make_array(mosaic_cols, mosaic_rows, /double)for ix = 0, mosaic_cols - 1 do xs[ix, *] = x_min + ix * raw_x_res + raw_x_res / 2.0dfor ix = 0, mosaic_rows - 1 do ys[*, ix] = y_max - ix * raw_y_res - raw_y_res / 2.0d; 注意: 非y_min + ix * raw_y_resproj_info = map_proj_init('Sinusoidal', sphere_radius=6371007.181d, /gctp, $ ; GCTP表示不使用IDL自带投影库进行而是基于GCTP进行center_longitude=0.0, false_easting=0.0, false_northing=0.0); 上述参数查看: https://modis-land.gsfc.nasa.gov/GCTP.htmllon_lat = map_proj_inverse(xs[valid_pixel_ix], ys[valid_pixel_ix], map_structure=proj_info)lon = lon_lat[0, *]lat = lon_lat[1, *]; GLT校正valid_pixel_cols_rows = array_indices(mosaic_imgs, valid_pixel_ix) ; 下标转换为行列号(shape=(2, None), 0为列, 1为行, None为有效像元个数); 输入的行列号数据集in_x = valid_pixel_cols_rows[0, *]in_y = valid_pixel_cols_rows[1, *]; 输出的行列号数据集out_x = floor((lon - (min(lon, /nan) - out_res / 2.0)) / out_res)out_y = floor(((max(lat, /nan) + out_res / 2.0) - lat) / out_res); 输出的行号和列号out_cols = ceil((max(lon, /nan) - min(lon, /nan)) / out_res)out_rows = ceil((max(lat, /nan) - min(lat, /nan)) / out_res); 取控制点子集(百分之~)pts_count = n_elements(valid_pixel_ix)sub_pts_count = floor(pts_count * sub_percent)sub_ix = floor(findgen(sub_pts_count) / sub_pts_count * pts_count)polywarp, in_x[sub_ix], in_y[sub_ix], out_x[sub_ix], out_y[sub_ix], degree, kx, ky ; 获取输入行列号与输出行列号的对应关系(系数)mosaic_imgs_warpped = poly_2d(mosaic_imgs, kx, ky, interp_code, out_cols, out_rows, missing=!values.F_NAN, $pixel_center=0.5) ; 基于对应关系进行输出; 地理结构体geo_info={$MODELPIXELSCALETAG: [out_res, out_res, 0.0], $ ; 分辨率MODELTIEPOINTTAG: [0.0, 0.0, 0.0, min(lon, /nan), max(lat, /nan), 0.0], $ ; 角点信息GTMODELTYPEGEOKEY: 2, $ ; 设置为地理坐标系GTRASTERTYPEGEOKEY: 1, $ ; 像素的表示类型, 北上图像(North-Up)GEOGRAPHICTYPEGEOKEY: 4326, $ ; 地理坐标系为WGS84GEOGCITATIONGEOKEY: 'GCS_WGS_1984', $GEOGANGULARUNITSGEOKEY: 9102} ; 单位为度; 输出out_path = out_dir + img_date + '.tiff'write_tiff, out_path, mosaic_imgs_warpped, geotiff=geo_info, /doubleprint, img_date, ' 时期镶嵌成功: ', systime(1) - cur_time, 's', format='%s%s%0.2f%s'endforeachprint, '总计用时: ', systime(1) - start_time, 's', format='%s%0.2f%s'
end
3.3 自定义时间类
时间类的定义相关概念和基本使用查看：https://blog.csdn.net/m0_63001937/article/details/133975751
本代码相关的date类：
function date::Init, _extra=ex; 初始化_ = self.idl_object::init() ; 初始化父类对象if isa(ex) then self.setproperty, _extra=exreturn, 1 ; 1表示成功实例化对象
endfunction date::strftime, format; 若不传入参数则默认该格式化字符串if ~keyword_set(format) then format='%Y-%M-%D'format = format.replace('%Y', string(self.year, format='%04d'))format = format.replace('%M', string(self.month, format='%02d'))format = format.replace('%D', string(self.day, format='%02d'))format = format.replace('%J', string(self.doy, format='%03d'))return, format
endfunction date::doy2ymd, year, doy, str=str, int=int; 更新if (keyword_set(year) and keyword_set(doy)) then beginself.year = yearself.doy = doyendif; 年积日转化为年月日old_date = imsl_datetodays(31, 12, self.year-1)imsl_daystodate, old_date + self.doy, day, month, yearself.day = dayself.month = monthself.year = yearreturn, self.to_ymd(str=str, int=int)
endfunction date::to_ymd, str=str, int=intif keyword_set(str) then return, string(self.year, self.month, self.day, format='%04d%02d%02d')if keyword_set(int) then return, long(string(self.year, self.month, self.day, format='%04d%02d%02d'))return, [self.year, self.month, self.day]
endfunction date::ymd2doy, year, month, day, str=str, int=int; 更新if (keyword_set(year) and keyword_set(month) and keyword_set(day)) then beginself.year = yearself.month = monthself.day = dayendif; 年月日转年积日now_days = imsl_datetodays(self.day, self.month, self.year)old_days = imsl_datetodays(31, 12, self.year - 1)self.doy = now_days - old_daysreturn, self.to_doy(str=str, int=int)
endfunction date::to_doy, str=str, int=intif keyword_set(str) then return, string(self.year, self.doy, format='%04d%03d')if keyword_set(int) then return, long(string(self.year, self.doy, format='%04d%03d'))return, [self.year, self.doy]
end; 更新年月日等属性
pro date::_update, ymd=ymd, doy=doyif keyword_set(ymd) then begin_ = self.ymd2doy(self.year, self.month, self.day)endifif keyword_set(doy) then begin_ = self.doy2ymd(self.year, self.doy)endif
end; 取值
pro date::GetProperty, year=year, month=month, day=day, doy=doy; 如果用户请求返回某一属性, 那么将其返回if arg_present(year) then year = self.yearIF arg_present(month) THEN month = self.monthIF arg_present(day) THEN day = self.dayif arg_present(doy) then doy = self.doy
end; 设置值
pro date::SetProperty, ymd=ymd, ydays=ydays; 如果用户传入属性, 那么就设置它if (isa(ymd, /integer) or isa(ymd, /string)) then beginif (isa(ymd, /integer)) then ymd = string(ymd, format='%08d')self.year = fix(strmid(ymd, 0, 4))self.month = fix(strmid(ymd, 4, 2))self.day = fix(strmid(ymd, 6, 2))self._update, /ymdendifif (isa(ydays, /integer) or isa(ydays, /string)) then beginif isa(ydays, /integer) then ydays = string(ydays, format='%07d')self.year = fix(strmid(ydays, 0, 4))self.doy = fix(strmid(ydays, 4, 3))self._update, /doyendif
end; 定义Date类
pro date__definestruct = {date, inherits idl_object, year: 0l, month: 0l, day: 0l, doy: 0l}
end
时间精力有限，代码各个细节虽很重要但无法一一说明。
这篇关于ENVI IDL：如何MODIS GRID产品进行批量镶嵌、重投影(GLT校正)？的文章就介绍到这儿，希望我们推荐的文章对编程师们有所帮助！
