遥感IDL语言保姆式提取GPM日降水量栅格数据的方法

数据源:NASA EARTHDATA 数据中心

地址：https://search.earthdata.nasa.gov/search?q=GPM&lat=36.99832109741731&long=-37.09165249454537&zoom=1

由于时间不够充裕所以没去做批处理只有单个的数据处理，大家可以自己做一下批处理，这样处理一堆批量的数据就快很多，如果有问题请给予指正~

空间分辨率:0.1度（约为10km*10km）

该数据优势：超强实时观测性（日降水最新数据距离现在时间前两天左右，半小时数据距离前一天左右) 观测范围广(一个数据观测全球)

数据挑选：由于GPM数据时间上分为每半小时/每天/每月，数据处理上分为早期/晚期/最终。看介绍由于多重数据的校正''最终''的精确度最高，但是最终截止于2021年9月还未更新，所以这里在数据选择上选择每日的晚期降水数据

2.查看数据集信息

不知道是不是因为这个数据是nc文件所以用hdf explorer 打不开，还是因为我的hdf explorer 本身有问题，不管了，用另外一种方法----ENVI  启动！

在这里也能较为清晰看到所包含的数据集

那我们应该怎样选择呢？

查阅官方文档....前为原版  后为翻译的中文版本

所以我们选择使用PrecipitationCal这个数据集

查阅属性--- （units--单位mm、long_name --长名称 Daily accumulated precipitation (combined microwave-IR) estimate 日累计降水量、_FillValue --填充值-9999....）

可以用得上的数据集为PrecipitationCal  lat  lon  即降水量数据 纬度、经度数据

提取步骤大白话来说无非就是：读降水数据-读经纬度数据-把经纬度拼成GLT再把GLT赋予降水数据-把赋予经纬度数据的降水数据进行重投影

这里要注意的是这个经纬度数据只是一个(1*1800 或者1*3600)的一维数组  简单来说就是行为1 列为1800或者3600的一条线，因此我们待会在处理数据的过程中先建立两个装经纬度数据的大数组--（行列号为1800*3600)。再把这两个代表为经纬度的数组输出为tiff和降水量数据输出的tiff做连接

输出的经度tiff（前）纬度tiff（后） 经度数据表示的是上面到下面出现数值的变化，纬度数据表示左右出现数值的变化，为什么和我们认知的经纬度数值变化正好相反？

看降水数据就能看明白了---因为降水数据本来也是倒过来的，当我们把经纬度数据赋予到降水量数据上，再进行重投影那么就会变得正常

废话少说直接上代码

即改即用：

代码：

pro get_pretem

  compile_opt idl2

  envi,/restore_base_save_files

  envi_batch_init

  file_name='E:\danzi\81\3B-DAY-L.MS.MRG.3IMERG.20230801-S000000-E235959.V06.nc4' 

  out_file='E:\danzi\81\'

  ncdf_id=ncdf_open(file_name);打开任意一个文件，并将文件保存成索引，因为所有数据的经纬度都一样

  ;处理经纬度信息获得像元尺度以及左上角像元经纬度

  att_name1='lon';我需要的变量名称

  att_name2='lat'

  att_name3='precipitationCal'

  varid1=ncdf_varid(ncdf_id,att_name1);根据文件索引，找到变量，并保存其变量索引

  varid2=ncdf_varid(ncdf_id,att_name2);根据文件索引，找到变量，并保存其变量索引

  varid3=ncdf_varid(ncdf_id,att_name3);根据文件索引，找到变量，并保存其变量索引

  ncdf_varget,ncdf_id,varid1,lon;获得变量 lon

  ncdf_varget,ncdf_id,varid2,lat;获得变量 lat

  ncdf_varget,ncdf_id,varid3,pre

  ncdf_attget,ncdf_id,varid3,'_FillValue',fv

;  lieshu=(size(lon))[1];列数

;  max_lon=max(lon);最大经度

;  min_lon=min(lon);最小经度

;  max_lat=max(lat);最大经度

;  min_lat=min(lat);最小经度

;  piexel_value=(max_lon-min_lon)/lieshu;像元值

;  ; 循环处理数据集

  pre=(pre ne fv[0])*pre

  lies=(size(lon))[1]

  hangs=(size(lat))[1]

  container1=fltarr(hangs,lies)

  container2=fltarr(hangs,lies)

  print,(size(container1))

  print,(size(lon))[1]

  for j=0,(size(container1))[1]-1 do begin

    for i=0,(size(container1))[2]-1 do begin

      container1[j,i]=lon[i]

    endfor

  endfor

  for j=0,(size(container2))[2]-1 do begin

    for i=0,(size(container2))[1]-1 do begin

      container2[i,j]=lat[i]

    endfor

  endfor

  out_name1=out_file+'lon.tiff'

  out_name2=out_file+'lat.tiff'

  out_name3=out_file+'pre.tiff'

  write_tiff,out_name1,container1,/float

  write_tiff,out_name2,container2,/float

  write_tiff,out_name3,pre,/float

  

  ;调用envi接口

  envi_open_file,out_name1,r_fid=x_fid

  envi_open_file,out_name2,r_fid=y_fid

  envi_open_file,out_name3,r_fid=data_fid

  ;输入glt栅格影像以及头文件输出路径

  out_name_glt=out_file+'preglt.img'

  out_name_glt_hdr=out_file+'preglt.hdr'

  ;定义输入输出坐标系

  input_proj=envi_proj_create(/geographic)

  output_proj=envi_proj_create(/geographic)

  envi_glt_doit,$

    x_fid=x_fid,y_fid=y_fid,x_pos=0,y_pos=0,i_proj=input_proj,$

    o_proj=output_proj,pixel_size=pixel_size,rotation=0.0,out_name=out_name_glt,r_fid=glt_id

  ;重投影

  out_name_geo=out_file+'pregeoref.img'

  out_name_geo_hdr=out_file+'pregeoref.hdr'

  envi_georef_from_glt_doit,$

    glt_fid=glt_id,pos=0,$

    fid=data_fid,$

    out_name=out_name_geo,background=0,r_fid=geo_id;赋予输出数据一个id号

  map_info=envi_get_map_info(fid=geo_id)

  geo_loc=map_info.(1)

  psize=map_info.(2)

  envi_file_query,geo_id,dims=data_dims

  target_data=envi_get_data(fid=geo_id,pos=0,dims=data_dims)

  geo_info={$

    MODELPIXELSCALETAG:[psize[0],psize[1],0.0],$

    MODELTIEPOINTTAG:[0.0,0.0,0.0,geo_loc[2],geo_loc[3],0.0],$

    GTMODELTYPEGEOKEY:2,$

    GTRASTERTYPEGEOKEY:1,$

    GEOGRAPHICTYPEGEOKEY:4326,$

    GEOGCITATIONGEOKEY:'GCS_WGS_1984',$

    GEOGANGULARUNITSGEOKEY:9102,$

    GEOGSEMIMAJORAXISGEOKEY:6378137.0,$

    GEOGINVFLATTENINGGEOKEY:298.257221}

  out_namedata=out_file+'81.tiff'

  write_tiff,out_namedata,target_data,/float,geotiff=geo_info

;  geo_info={$

;    MODELPIXELSCALETAG:[piexel_value,piexel_value,0.0],$

;    MODELTIEPOINTTAG:[0.0,0.0,0.0,min_lon,max_lat,0.0],$

;    GTMODELTYPEGEOKEY:2,$

;    GTRASTERTYPEGEOKEY:1,$

;    GEOGRAPHICTYPEGEOKEY:4326,$

;    GEOGCITATIONGEOKEY:'GCS_WGS_1984',$

;    GEOGANGULARUNITSGEOKEY:9102,$

;    GEOGSEMIMAJORAXISGEOKEY:6378137.0,$

;    GEOGINVFLATTENINGGEOKEY:298.257221}

;  out_name1=out_file+'pre.tiff'

;  write_tiff,out_name1,pre,/float,geotiff=geo_info

envi_file_mng,id=x_fid,/remove

envi_file_mng,id=y_fid,/remove

envi_file_mng,id=data_fid,/remove

envi_file_mng,id=glt_id,/remove

envi_file_mng,id=geo_id,/remove

file_delete,[out_name1,out_name2,out_name3,out_name_geo,out_name_geo_hdr,out_name_glt,out_name_glt_hdr]

envi_batch_exit,/no_confirm

end

拿2023年8月1日 华北地区出现洪涝灾害这一天数据为例效果如下：