IDL常规数据的读取

空格 或 制表符 分隔数据读取方法

1.打开文件(OPENR)
2.获取文件行数(FILE LINES)
3.获取文件列数(手动查看或编程获取)
4.建立对应行列数的数组
5.读取数据(READF)、（SKIP_LUN）
6.释放设备单元（FREE_LUN）

为什么有自动确定行的函数，却没有自动确定列的函数呢？
因为行数可以通过 换行符 进行判断，但是列，可能的分割方式有很多，需要自己去确定。

可以看到可以用空格做分割，读取第一行后，用空格做分割后，有几部分，就有几列

pro txt_read
  file_name = 'E:\Data\IDL\chapter_4\2013_year_aop.txt'
  openr,1,file_name
  str = ' ';定义一个空的字符串,因为只有一行，所以只能读一行
  readf,1,str;定义的变量有几行，readf就只能读几行
  print,str
;  str = strarr(1,2);定义一个2行的字符串
;  readf,1,str;定义的变量有2行，readf就只能读2行
  ;读取之后，列之间是以空格为分割，所以以下才做能分割出多少个元素，就是有几列
  parameters = n_elements(strsplit(str,' '))  
  parameters = n_elements(strsplit(str,' '))
  print,parameters
  print,strsplit(str,' ',/extract)
  line = file_lines(file_name)
  
  ;print,str
  line = file_lines(file_name)
  free_lun,1
end

这样就可以创建需要的数组大小，进行储存了。

  line = file_lines(file_name) - 1 ;是因为表头（第一行）是不要的
  data = fltarr(parameters,line)
  readf,1,data

CSV数据读取，逗号 分隔数据读取方法

  file_name = 'E:\Data\IDL\chapter_4\air_quality.csv'
  data = read_csv(file_name,header = par_name)  ;默认不读表头，如果想要第一行表头，则需要加关键字,header = 
  print,data ; 按列打出
  print,par_name ; 表头
  help,data ; 此时data是一个结构体

此时data是结构体，我们可以靠表头得知我们需要的内容对应的是哪一个内容

  data = read_csv(file_name,header = par_name)  ;默认不读表头，如果想要第一行表头，则需要加关键字,header = 
  print,data ; 按列打出
  print,par_name ; 表头
  help,data ; 此时data是一个结构体
  
  pos = where(par_name eq 'SO2')
  print,pos
  so2 = data.(pos)
  print,so2

文本数据可视化（txt转成geotiff）

核心：
数据点的空间位置确定
问题分解：
1.图幅大小(与经纬度范围、像元分辨率相关)·
2.空间位置
3.计算与填入
4.空值填补

当一个函数有多个返回值的时候，可以写成pro函数。

pro ..........
......
end

pro后面不仅仅有输入参数，也要有输出参数。会自动识别哪些是输入参数，哪些是输出参数。

自动重投影函数：

pro geo_box_generating,longitude,latitude,target_data,resolution,geo_box,geo_box_out,geo_info
  lon_min=min(longitude)
  lon_max=max(longitude)
  lat_min=min(latitude)
  lat_max=max(latitude)
  
  geo_box_col=ceil((lon_max-lon_min)/resolution)
  geo_box_line=ceil((lat_max-lat_min)/resolution)
  geo_box=fltarr(geo_box_col,geo_box_line)
  
  geo_box_col_pos=floor((longitude-lon_min)/resolution)
  geo_box_line_pos=floor((lat_max-latitude)/resolution)
  geo_box[geo_box_col_pos,geo_box_line_pos]=target_data
  
;此时完成了经纬度重投影，还需要填补空缺值，并将图片第一行、最后一行、第一列、最后一列保存。

;补全四周行列
  
  geo_box_copy=fltarr(geo_box_col+2,geo_box_line+2)
  geo_box_copy[1:geo_box_col,1:geo_box_line]=geo_box
  
;补全空缺值
  geo_box_out=fltarr(geo_box_col,geo_box_line)
  for geo_box_col_i=1,geo_box_col do begin
    for geo_box_line_i=1,geo_box_line do begin
      if geo_box_copy[geo_box_col_i,geo_box_line_i] eq 0.0 then begin
        temp_window=geo_box_copy[geo_box_col_i-1:geo_box_col_i+1,geo_box_line_i-1:geo_box_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
          geo_box_out[geo_box_col_i-1,geo_box_line_i-1]=temp_window_sum/temp_window_num
        endif
      endif else begin
        geo_box_out[geo_box_col_i-1,geo_box_line_i-1]=geo_box_copy[geo_box_col_i,geo_box_line_i]
      endelse
    endfor
  endfor
  
  geo_info={$
      MODELPIXELSCALETAG:[resolution,resolution,0.0],$
      MODELTIEPOINTTAG:[0.0,0.0,0.0,lon_min,lat_max,0.0],$
      GTMODELTYPEGEOKEY:2,$
      GTRASTERTYPEGEOKEY:1,$
      GEOGRAPHICTYPEGEOKEY:4326,$
      GEOGCITATIONGEOKEY:'GCS_WGS_1984',$
      GEOGANGULARUNITSGEOKEY:9102,$
      GEOGSEMIMAJORAXISGEOKEY:6378137.0,$
      GEOGINVFLATTENINGGEOKEY:298.25722}
end

此时，主函数为：

pro txt_read
  filename='E:/Data/IDL/chapter_4/2013_year_aop.txt'
  openr,1,filename
  str=' '
  readf,1,str
  parameters=strsplit(str,' ',/extract)
  par_n=n_elements(parameters)
  line=file_lines(filename)-1
  data=fltarr(par_n,line)
  readf,1,data
  free_lun,1
  
  longitude=data[0,*]
  latitude=data[1,*]
  resolution=0.18
  
  for par_i=2,par_n-1 do begin
    out_tiff='E:/Data/IDL/chapter_4/out_tiff/'+parameters[par_i]+'.tiff'
    target_data=data[par_i,*]
    geo_box_generating,longitude,latitude,target_data,resolution,geo_box,geo_box_out,geo_info
    write_tiff,out_tiff,geo_box_out,/float,geotiff=geo_info
  endfor
end

为什么有多个输出？
是方便我们去调用不同的输出去做后续应用。

数据储存为hdf5文件

有时候要求不要 tiff 而是需要输出 hdf 文件。


1.这里的 *.tiff的含义是 所有的 tiff 文件

  file_list = file_search(input_directory,'*.tiff') 

2.这里的month *.tiff的含义是 所有month开头的 tiff 文件。

  file_list = file_search(input_directory,'month*.tiff') 

3.这里的 * month *.tiff 含义是所有命名当中有month的文件，开头结尾不管，只要中间有month就读取。

  file_list = file_search(input_directory,'*month*.tiff') 

hdf5的创建与关闭

  file_id = h5f_create(hdf5_out_name)
  
  h5f_close,file_id

创建一个通用属性信息：

  file_id = h5f_create(hdf5_out_name)
    att_value_temp = 'OMI monthly NO2 product of 2017~2018' ; 写一个属性
    datatype_id = h5t_idl_create(att_value_temp) ;虽然都是字符串，但hdf5和计算机底层的字符串类型是不一样的，所以需要一个转换
    dataspace_id = h5s_create_simple([1]);确定需要多大的空间去存储，[1]代表用一个元素的数组去存储
    attr_id = h5a_create(file_id,'Description',datatype_id,dataspace_id);datatype_id数据类型，dataspace_id需要占多大空间
    h5a_write,attr_id,att_value_temp
   
  h5f_close,file_id

创建分组数据，在组下面加属性

    data_group_id = h5g_create(file_id,'OMI Monthly Data')

file_id = h5f_create(hdf5_out_name)
    att_value_temp = 'OMI monthly NO2 product of 2017~2018' ; 写一个属性
    datatype_id = h5t_idl_create(att_value_temp) ;虽然都是字符串，但hdf5和计算机底层的字符串类型是不一样的，所以需要一个转换
    dataspace_id = h5s_create_simple([1]);确定需要多大的空间去存储，[1]代表用一个元素的数组去存储
    attr_id = h5a_create(file_id,'Description',datatype_id,dataspace_id);datatype_id数据类型，dataspace_id需要占多大空间
    ;h5a_create是创建数据属性
    
    h5a_write,attr_id,att_value_temp
    
    data_group_id = h5g_create(file_id,'OMI Monthly Data')
    for file_i = 0,n_elements(file_list)-1 do begin
      target_data = read_tiff(file_list[file_i],geotiff = geo_info)
      data_size = size(target_data)
      data_col = data_size[1]
      data_line = data_size[2]
      resorlution_tag = geo_info.(0)
      geo_tag = geo_info.(1)
      dataset_name = file_basename(file_list[file_i],'.tiff');以文件名当做数据集名字，并把后缀去除
      
      data_type_id = h5t_idl_create(target_data)
      dataspace_id =h5s_create_simple([data_col,data_line]) 
      dataset_id = h5d_create(data_group_id,dataset_name,data_type_id,dataspace_id)
      ;h5d_create创建数据集
      h5d_write,dataset_id,target_data
      
    endfor

数据压缩

此时，没有数据压缩，如果想使得数据压缩，那么可以加2个东西：

      dataset_id = h5d_create(data_group_id,dataset_name,data_type_id,dataspace_id,$
        chunk_dimension=[data_col,data_line],GZIP = 9)

1.h5d_create创建数据集
2.chunk_dimension=[data_col,data_line]代表被压缩数据的大小，一般是全部压缩
3.GZIP压缩等级，9是压缩最高等级

在数据集下面加属性

向往那个下面写，就往哪个id下面创建。
注意：同一个对象同样的属性只能被写入一次。
部分代码：

      ;在数据集下面加属性
      att_value_temp = 0 
      datatype_id = h5t_idl_create(att_value_temp) ;虽然都是字符串，但hdf5和计算机底层的字符串类型是不一样的，所以需要一个转换
      dataspace_id = h5s_create_simple([1]);确定需要多大的空间去存储，[1]代表用一个元素的数组去存储
      attr_id = h5a_create(dataset_id,'_FillValue',datatype_id,dataspace_id);datatype_id数据类型，dataspace_id需要占多大空间
      h5a_write,attr_id,att_value_temp

全部代码：

pro hdf5_writing
  input_directory = 'E:/Data/IDL/chapter_2/NO2/average/'
  output_directory = 'E:/Data/IDL/chapter_2/NO2/average/hdf5/'
  dir_test = file_test(output_directory,/directory)
  if dir_test eq 0 then file_mkdir,output_directory ;检测文件是否存在，若不存在则创建
  hdf5_out_name = output_directory + 'OMI_NO2_month_product.h5'
  file_list = file_search(input_directory,'*month*.tiff') 
  
  file_id = h5f_create(hdf5_out_name)
    att_value_temp = 'OMI monthly NO2 product of 2017~2018' ; 写一个属性
    datatype_id = h5t_idl_create(att_value_temp) ;虽然都是字符串，但hdf5和计算机底层的字符串类型是不一样的，所以需要一个转换
    dataspace_id = h5s_create_simple([1]);确定需要多大的空间去存储，[1]代表用一个元素的数组去存储
    attr_id = h5a_create(file_id,'Description',datatype_id,dataspace_id);datatype_id数据类型，dataspace_id需要占多大空间
    ;h5a_create是创建数据属性
    
    h5a_write,attr_id,att_value_temp
    
    data_group_id = h5g_create(file_id,'OMI Monthly Data')
    for file_i = 0,n_elements(file_list)-1 do begin
      target_data = read_tiff(file_list[file_i],geotiff = geo_info)
      data_size = size(target_data)
      data_col = data_size[1]
      data_line = data_size[2]
      resorlution_tag = geo_info.(0)
      geo_tag = geo_info.(1)
      dataset_name = file_basename(file_list[file_i],'.tiff');以文件名当做数据集名字，并把后缀去除
      
      data_type_id = h5t_idl_create(target_data)
      dataspace_id =h5s_create_simple([data_col,data_line]) 
      dataset_id = h5d_create(data_group_id,dataset_name,data_type_id,dataspace_id,$
        chunk_dimension=[data_col,data_line],GZIP = 9)
      h5d_write,dataset_id,target_data
      ;h5d_create创建数据集
      ;chunk_dimension=[data_col,data_line]代表被压缩数据的大小，一般是全部压缩
      ;GZIP压缩等级，9是压缩最高等级
      
      ;在数据集下面加属性
      att_value_temp = 0 
      datatype_id = h5t_idl_create(att_value_temp) ;虽然都是字符串，但hdf5和计算机底层的字符串类型是不一样的，所以需要一个转换
      dataspace_id = h5s_create_simple([1]);确定需要多大的空间去存储，[1]代表用一个元素的数组去存储
      attr_id = h5a_create(dataset_id,'_FillValue',datatype_id,dataspace_id);datatype_id数据类型，dataspace_id需要占多大空间
      h5a_write,attr_id,att_value_temp    
    endfor 
  h5f_close,file_id
 
end

创建hdf的全部代码：

pro tiff_to_hdf5
  input_directory='O:/coarse_data/chapter_2/NO2/average/'
  output_directory='O:/coarse_data/chapter_2/NO2/average/hdf5/'
  dir_test=file_test(output_directory,/directory)
  if dir_test eq 0 then begin
    file_mkdir,output_directory
  endif
  hdf5_out_name=output_directory+'OMI_NO2_Month_Product.h5'
  file_list=file_search(input_directory,'month*.tiff')
  
  file_id=h5f_create(hdf5_out_name)
  att_value_temp='OMI monthly NO2 product of 2017 - 2018'
  datatype_id=h5t_idl_create(att_value_temp)
  dataspace_id=h5s_create_simple([1])
  attr_id=h5a_create(file_id,'Description',datatype_id,dataspace_id)
  h5a_write,attr_id,att_value_temp
  
  data_group_id=h5g_create(file_id,'OMI monthly data')
  
  for file_i=0,n_elements(file_list)-1 do begin
    target_data=read_tiff(file_list[file_i],geotiff=geo_info)
    data_size=size(target_data)
    data_box_col=data_size[1]
    data_box_line=data_size[2]
    resolution_tag=geo_info.(0)
    geo_tag=geo_info.(1)    
    dataset_name=file_basename(file_list[file_i],'.tiff')
    
    datatype_id=h5t_idl_create(target_data)
    dataspace_id=h5s_create_simple([data_box_col,data_box_line])
    dataset_id=h5d_create(data_group_id,dataset_name,datatype_id,dataspace_id,$
      chunk_dimensions=[data_box_col,data_box_line],gzip=9)
    h5d_write,dataset_id,target_data
    
    att_value_temp=0
    datatype_id=h5t_idl_create(att_value_temp)
    dataspace_id=h5s_create_simple([1])
    attr_id=h5a_create(dataset_id,'Fill_Value',datatype_id,dataspace_id)
    h5a_write,attr_id,att_value_temp

    att_value_temp=1
    datatype_id=h5t_idl_create(att_value_temp)
    dataspace_id=h5s_create_simple([1])
    attr_id=h5a_create(dataset_id,'Scale_Factor',datatype_id,dataspace_id)
    h5a_write,attr_id,att_value_temp   
    
    att_value_temp='mol / km2'
    datatype_id=h5t_idl_create(att_value_temp)
    dataspace_id=h5s_create_simple([1])
    attr_id=h5a_create(dataset_id,'Unit',datatype_id,dataspace_id)
    h5a_write,attr_id,att_value_temp
  endfor
  ;geo_group_id=h5g_create(data_group_id,'geo data')
  geo_group_id=h5g_create(file_id,'geo data')
  datatype_id=h5t_idl_create(resolution_tag[0:1])
  dataspace_id=h5s_create_simple([2])
  dataset_id=h5d_create(geo_group_id,'grid_space',datatype_id,dataspace_id)
  h5d_write,dataset_id,resolution_tag[0:1]
  
  datatype_id=h5t_idl_create(geo_tag[3:4])
  dataspace_id=h5s_create_simple([2])
  dataset_id=h5d_create(geo_group_id,'ul_corner',datatype_id,dataspace_id)
  h5d_write,dataset_id,geo_tag[3:4]
  
  h5a_close,attr_id
  h5t_close,datatype_id
  h5s_close,dataspace_id
  h5d_close,dataset_id
  h5g_close,gid
  h5f_close,file_id
end

以上还对geo_info结构体所需要的信息进行了一个储存。

ENVI数据如何去读取

envi数据全部有用信息全部在hdr文件下面。

interleave是envi遥感图像的排列方式
byte order是大小端的概念
map info是地理坐标信息
data type 数据类型,1 就是byt型
建立的时候需要，行列号以及类型，和目标数据一致，必须一致才能读对。

image(data[*,*,0]);读取第一层的数据

但是本身image读取可能会存在一定的颠倒，所以需一个旋转。

image(rotate(data[*,*,0],7))