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.0
    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 ; 对应行列号上应放置的像素是哪一个

原理参考如下：



填充条件：
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 begin
      for data_box_geo_line_i=1,data_box_geo_line-2 do begin
        if data_box_geo[data_box_geo_col_i,data_box_geo_line_i] eq -9999.0 then begin
          temp_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_window
          temp_window_sum=total(temp_window)
          temp_window_num=total(temp_window gt 0.0)
          if (temp_window_num gt 3) then begin
            data_box_geo_out[data_box_geo_col_i,data_box_geo_line_i]=temp_window_sum/temp_window_num
          endif
        endif else begin
          data_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]
        endelse
      endfor
    endfor

对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_geographic
  input_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 begin
    file_mkdir,output_directory
  endif
  file_list = file_search(input_directory,'*.hdf')
  file_n = n_elements(file_list)
  ;print, file_n
  
  for file_i = 0, file_n - 1 do begin ;file_n - 1 而不是file_n
    start_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, data
    
    index = hdf_sd_attrfind(sds_id, 'scale_factor') ; 为什么这次没有计算——fillvalue? 
    hdf_sd_attrinfo,sds_id, index, COUNT=cali_num,DATA=cali_scale
    data = data * cali_scale[0]
    hdf_sd_endaccess,sds_id
    hdf_sd_end,sd_id
    
    data_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 begin
      proj_x[col_i,*] = ul_prj_x + (sin_resolution * col_i) ; 利用 左上坐标、像元分辨率，计算每个点坐标
    endfor
    
    for line_i = 0, data_size[2] - 1 do begin
      proj_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_pos
    geo_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 begin
      for data_box_geo_line_i=1,data_box_geo_line-2 do begin
        if data_box_geo[data_box_geo_col_i,data_box_geo_line_i] eq -9999.0 then begin
          temp_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_window
          temp_window_sum=total(temp_window)
          temp_window_num=total(temp_window gt 0.0)
          if (temp_window_num gt 3) then begin
            data_box_geo_out[data_box_geo_col_i,data_box_geo_line_i]=temp_window_sum/temp_window_num
          endif
        endif else begin
          data_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]
        endelse
      endfor
    endfor
    
    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}
    
    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

原文地址：https://blog.csdn.net/Aigcl/article/details/144300222

免责声明：本站文章内容转载自网络资源，如本站内容侵犯了原著者的合法权益，可联系本站删除。更多内容请关注自学内容网（zxcms.com）！