IDL学习笔记（五）MODIS数据（Grid）

正弦投影 是以 米 为单位的
经纬度网格 是以 度 为单位的

但是转换之后，不会一一对应，所以需要对中间空缺位置需要进行一个填补。

核心问题:

把一个点从一个空间参考系放到另一个空间参考系

对应的编程关键:

1.数组A的元素应该放到数组B的哪个位置?

2.数组B的行列数量?

3.数组B的空值如何填补？

MODIS Grid数据的重投影

1.对应的专业问题:

数组A的元素应该放到数组B的哪个位置投影

2.坐标与经纬度坐标的转换：

数组B的行列数量:经纬度格网的范围

3.数组B的空值如何填补:

移动窗口均值计算



首先，确定，经纬度最值。（数组行列数:原图像包含的经纬度范围最值 ； 重投影结果分辨率设置）

重投影结果分辨率设置，不能比原来更细，因为这样的话，空缺会更多，将带来更大误差。

每个像元坐标从哪里来？

知道第一列坐标，那么第二列坐标 就是 第一列坐标 + 像元分辨率。所以其余坐标就可以利用像元分辨率进行 加减乘除 运算 得知。

像元分辨率没有提到，需要自己计算。利用左上角、右下角坐标，以及行数、列数，进行计算得到。

再找到，我们需要的“一行”信息，再定位到（）取出里面的坐标信息。为什么不用strmid？因为（）内的信息可能是不固定的，长度不唯一。

根据前面的关键字，通过strpos定位到起始位置和结束位置，获取总共子串长度，然后就可以用strsplit分割函数获取所需数值。

strsplit(ul_info, '=(,)', /extract)

1.'=(,)'只要遇到引号里的任何一个内容就把他分开
2./extract 如果需要提取内容，则输出的是分割的下标。所以需要提取信息的话，需要添加关键字。这时候，提取结果是一个 字符型 数组，里面包含坐标信息。

    sd_id = hdf_sd_start(file_list[file_i],/read)gindex = hdf_sd_attrfind(sd_id,'StructMetadata.0')hdf_sd_attrinfo, sd_id, gindex, data=metadata;左上角坐标获取ul_start_pos = strpos(metadata,'UpperLeftPointMtrs')ul_end_pos = strpos(metadata,'LowerRightMtrs')ul_info = strmid(metadata,ul_start_pos, ul_end_pos-ul_start_pos)ul_info_spl = strsplit(ul_info,'=(,)',/extract) ; 加 /extract 才返回的是内容，而不是切割点的索引ul_prj_x = double(ul_info_spl[1])ul_prj_y = double(ul_info_spl[2]);右下角坐标获取lr_start_pos = strpos(metadata,'LowerRightMtrs')lr_end_pos = strpos(metadata,'Projection')lr_info = strmid(metadata,lr_start_pos,le_end_pos - le_start_pos)lr_split = strsplit(lr_info,'=(,)',/extract) ;这里一开始忘添加关键字了，导致后面查报错查了很久lr_prj_x = double(lr_split[1])lr_prj_y = double(lr_split[2])

	;map_proj_inverse 用于将 投影坐标 转为 经纬度坐标sin_prj = map_proj_init('sinusoidal', /gctp, sphere_radius = 6371007.181,center_longitude=0.0,false_easting=0.0,false_northing=0.0);投影坐标参数，从数据里得知，geo_loc = map_proj_inverse(proj_x, proj_y, map_strcture = sin_prj); 进行坐标转换，需要知道原来的投影的坐标

 	lon_min = min(geo_x)lon_max = max(geo_x)lat_min = min(geo_y)lat_max = max(geo_y)geo_resolution = 0.01 ;1公里的空间分辨率，对应的大概是0.01度data_box_geo_col = ceil((lon_max - lon_min)/geo_resolution) ; ceil 防止无法整除带来的误差，小数部分需要向上取整，保证数据全部存放，没有丢失部分data_box_geo_line = ceil((lat_max - lat_min)/geo_resolution)data_box_geo = fltarr(data_box_geo_col, data_box_geo_line) ; 初始化data_box_geo[*,*] = -9999.0data_box_geo_col_pos = floor((geo_x - lon_min)/geo_resolution) ; 归到哪一列，需要向下取整，因为data_box_geo_line_pos = floor((geo_y - lat_min)/geo_resolution)

对应位置存放：

data_box_geo[data_box_geo_col_pos, data_box_geo_line_pos] = data ; 对应行列号上应放置的像素是哪一个

原理参考如下：



填充条件：
1.同时遇到_FillValue（这里是0）和-9999（初始值），要忽略，不参与运算。
2.有N个以上有效值（例如3）的时候才进行填充，否则就是不填充.不填充结果也设置为0
3.不需要填充的数组，保持不动

对空缺值进行填充（平滑）：

data_box_geo_out = fltarr(data_box_geo_col, data_box_geo_line)for data_box_geo_col_i=1,data_box_geo_col-2 do beginfor data_box_geo_line_i=1,data_box_geo_line-2 do beginif data_box_geo[data_box_geo_col_i,data_box_geo_line_i] eq -9999.0 then begintemp_window=data_box_geo[data_box_geo_col_i-1:data_box_geo_col_i+1,data_box_geo_line_i-1:data_box_geo_line_i+1]temp_window=(temp_window gt 0.0)*temp_windowtemp_window_sum=total(temp_window)temp_window_num=total(temp_window gt 0.0)if (temp_window_num gt 3) then begindata_box_geo_out[data_box_geo_col_i,data_box_geo_line_i]=temp_window_sum/temp_window_numendifendif else begindata_box_geo_out[data_box_geo_col_i,data_box_geo_line_i]=data_box_geo[data_box_geo_col_i,data_box_geo_line_i]endelseendforendfor

对geo结构体进行定义

    geo_info = {$MODELPIXELSCALETAG:[geo_resolution,geo_resolution,0.0],$;x、y、z方向的像元分辨率MODELTIEPOINTTAG:[0.0,0.0,0.0,lon_min,lat_max,0.0],$;坐标转换信息，前三个0.0代表栅格图像上的第0,0,0个像元位置（z方向一般不存在），后面-180.0代表x方向第0个位置对应的经度是-180.0度，90.0代表y方向第0个位置对应的纬度是90.0度GTMODELTYPEGEOKEY:2,$GTRASTERTYPEGEOKEY:1,$GEOGRAPHICTYPEGEOKEY:4326,$GEOGCITATIONGEOKEY:'GCS_WGS_1984',$GEOGANGULARUNITSGEOKEY:9102,$GEOGSEMIMAJORAXISGEOKEY:6378137.0,$GEOGINVFLATTENINGGEOKEY:298.25722}

输出一个文件后，打印输出的提示信息，可以都设置为以下格式：

     print,'Reprojection time consuming of file '+file_basename(file_list[file_i])+':'+strcompress(string(end_time-start_time))+' s.'

完整代码：

pro modis_sinusoidal_to_geographicinput_directory = 'E:\Data\IDL\chapter_3\modis_grid\'output_directory = 'E:\Data\IDL\chapter_3\modis_grid\geo_out\'directory_exist = file_test(output_directory,/directory)if directory_exist eq 0 then beginfile_mkdir,output_directoryendiffile_list = file_search(input_directory,'*.hdf')file_n = n_elements(file_list);print, file_nfor file_i = 0, file_n - 1 do begin ;file_n - 1 而不是file_nstart_time = systime(1)result_name = output_directory + file_basename(file_list[file_i],'.hdf') + '_geo.tiff'sd_id = hdf_sd_start(file_list[file_i],/read)gindex = hdf_sd_attrfind(sd_id,'StructMetadata.0')hdf_sd_attrinfo, sd_id, gindex, data=metadata;左上角坐标获取ul_start_pos = strpos(metadata,'UpperLeftPointMtrs')ul_end_pos = strpos(metadata,'LowerRightMtrs')ul_info = strmid(metadata,ul_start_pos, ul_end_pos-ul_start_pos)ul_info_spl = strsplit(ul_info,'=(,)',/extract) ; 加 /extract 才返回的是内容，而不是切割点的索引ul_prj_x = double(ul_info_spl[1])ul_prj_y = double(ul_info_spl[2]);右下角坐标获取lr_start_pos = strpos(metadata,'LowerRightMtrs')lr_end_pos = strpos(metadata,'Projection')lr_info = strmid(metadata,lr_start_pos,lr_end_pos - lr_start_pos)lr_info_spl=strsplit(lr_info,'=(,)',/extract)lr_prj_x = double(lr_split[1])lr_prj_y = double(lr_split[2])sds_index = hdf_sd_nametoindex(sd_id,'LST_Day_1km')sds_id = hdf_sd_select(sd_id,sds_index)hdf_sd_getdata, sds_id, dataindex = hdf_sd_attrfind(sds_id, 'scale_factor') ; 为什么这次没有计算——fillvalue? hdf_sd_attrinfo,sds_id, index, COUNT=cali_num,DATA=cali_scaledata = data * cali_scale[0]hdf_sd_endaccess,sds_idhdf_sd_end,sd_iddata_size = size(data)sin_resolution = (lr_prj_x - ul_prj_x)/(data_size[1]) ; 计算原始正弦投影的 像元分辨率 proj_x = fltarr(data_size[1], data_size[2]) ; 初始化数组，存放坐标proj_y = fltarr(data_size[1], data_size[2])for col_i = 0, data_size[1] -1 do beginproj_x[col_i,*] = ul_prj_x + (sin_resolution * col_i) ; 利用 左上坐标、像元分辨率，计算每个点坐标endforfor line_i = 0, data_size[2] - 1 do beginproj_y[*, line_i] = ul_prj_y - (sin_resolution * line_i)endfor;map_proj_inverse 用于将 投影坐标 转为 经纬度坐标sin_prj = map_proj_init('sinusoidal', /gctp, sphere_radius = 6371007.181,center_longitude=0.0,false_easting=0.0,false_northing=0.0);投影坐标参数，从数据里得知，geo_loc = map_proj_inverse(proj_x, proj_y, map_structure = sin_prj); 进行坐标转换，需要知道原来的投影的坐标;geo是2列1440000行的数组，第一列是精度，第二列是纬度geo_x = geo_loc[0, *] ; x方向 对应应放入  data_box_geo_col_posgeo_y = geo_loc[1, *] ; y方向 对应应放入  data_box_geo_line_pos;到此为止，得知了原图上经纬度坐标。lon_min = min(geo_x)lon_max = max(geo_x)lat_min = min(geo_y)lat_max = max(geo_y)geo_resolution = 0.01 ;1公里的空间分辨率，对应的大概是0.01度data_box_geo_col = ceil((lon_max - lon_min)/geo_resolution) ; ceil 防止无法整除带来的误差，小数部分需要向上取整，保证数据全部存放，没有丢失部分data_box_geo_line = ceil((lat_max - lat_min)/geo_resolution)data_box_geo = fltarr(data_box_geo_col, data_box_geo_line) ; 初始化data_box_geo[*,*] = -9999.0 ; 为什么赋值-9999.0 ？是为了区分是原本数据缺失值 还是 本身就有的0值（_FillValue）,判断是否要做填充data_box_geo_col_pos = floor((geo_x - lon_min)/geo_resolution) ; 归到哪一列，需要向下取整，因为data_box_geo_line_pos = floor((geo_y - lat_min)/geo_resolution)data_box_geo[data_box_geo_col_pos, data_box_geo_line_pos] = data ; 对应行列号上应放置的像素是哪一个;对空缺值做一个填补操作----平滑data_box_geo_out = fltarr(data_box_geo_col, data_box_geo_line)for data_box_geo_col_i=1,data_box_geo_col-2 do beginfor data_box_geo_line_i=1,data_box_geo_line-2 do beginif data_box_geo[data_box_geo_col_i,data_box_geo_line_i] eq -9999.0 then begintemp_window=data_box_geo[data_box_geo_col_i-1:data_box_geo_col_i+1,data_box_geo_line_i-1:data_box_geo_line_i+1]temp_window=(temp_window gt 0.0)*temp_windowtemp_window_sum=total(temp_window)temp_window_num=total(temp_window gt 0.0)if (temp_window_num gt 3) then begindata_box_geo_out[data_box_geo_col_i,data_box_geo_line_i]=temp_window_sum/temp_window_numendifendif else begindata_box_geo_out[data_box_geo_col_i,data_box_geo_line_i]=data_box_geo[data_box_geo_col_i,data_box_geo_line_i]endelseendforendforgeo_info = {$MODELPIXELSCALETAG:[geo_resolution,geo_resolution,0.0],$;x、y、z方向的像元分辨率MODELTIEPOINTTAG:[0.0,0.0,0.0,lon_min,lat_max,0.0],$;坐标转换信息，前三个0.0代表栅格图像上的第0,0,0个像元位置（z方向一般不存在），后面-180.0代表x方向第0个位置对应的经度是-180.0度，90.0代表y方向第0个位置对应的纬度是90.0度GTMODELTYPEGEOKEY:2,$GTRASTERTYPEGEOKEY:1,$GEOGRAPHICTYPEGEOKEY:4326,$GEOGCITATIONGEOKEY:'GCS_WGS_1984',$GEOGANGULARUNITSGEOKEY:9102,$GEOGSEMIMAJORAXISGEOKEY:6378137.0,$GEOGINVFLATTENINGGEOKEY:298.25722}write_tiff, result_name, data_box_geo_out,/float, geo_tiff = geo_info end_time = systime(1)print,'Reprojection time consuming of file '+file_basename(file_list[file_i])+':'+strcompress(string(end_time-start_time))+' s.'endfor 
end