Showing posts with label IDL. Show all posts
Showing posts with label IDL. Show all posts

Wednesday, May 10, 2017

IDL: 关闭所有文件

代码示例： 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
;; get all files IDs
fids = envi_get_file_ids()
size = size(fids)
length = size[1]

for j = 0L, length-1, 1 do begin
      ;; to delete the specified file from the hard disk
      ENVI_FILE_MNG, id = fids[j], /DELETE
      ;; to remove the specified file from within ENVI Classic
      ;;ENVI_FILE_MNG, id = fids[j], /REMOVE
endfor
Posted by at No comments:
Labels: ,

Monday, June 6, 2016

IDL: MVC

Introduction

MVC（Maximum-value composite）方法计算得出多个矩阵之间对应位置的最大值，空间数据处理中比较常用，尤其是NDVI。这里介绍两个版本的IDL语言MVC函数，其区别在于计算过程的内存占用情况，若输入N个文件，第一版内存占用为N，第二版则仅有1/N，后者相比前者对内存要求较小，适用于单个文件较大的情况。

第一版

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
;;Created by LI Xu
;;20 November, 2013
;;Version 1.0

;;Discrption
;;MVC, the maximum-value composite method
;;infullname: string array, consist of N files with full path
;;strSaveDir: string, the dirctory of result
;;strflag: string, flag to month or year , etc..


;+
; :Description:
;    Describe the procedure.
;
; :Params:
;    infullname
;    strSaveDir
;    strflag
;
;
;
; :Author: LiXu
;-
pro MVC, infullname, strSavepath

  
  ;;Read one by one
  ok=Query_Tiff(infullname[0], info)
  Cols=info.dimensions[0]
  Rows=info.dimensions[1]
  N_Dim=size(infullname, /n_elements)
  Pics=DblArr(N_Dim, Cols, Rows)
  for ii=0, N_Dim-1 do Begin
    Pics[ii, *, *]=Read_Tiff(infullname[ii], GEOTIFF=GEOVAR)   
  endfor
  
  ;;Run > 
  Result=Pics[0, *, *]
  ;;Loop >
  for jj=1, N_Dim-1 do Begin
    Result=Result>Pics[jj, *, *]   
  endfor
  
  ;;Write a result tiff
  Write_Tiff, strSavePath, Result, GEOTIFF=GEOVAR, /float
  print, strSavePath+' done!'
end

第二版

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
;; Created by LI Xu
;; Version 1.0
;; June 6, 2016

;; Description:
;; MVC Method

;; If you have any question about this code,
;; please do not hesitate to contact me via E-mail:
;; jeremy456@163.com

;; Blog:
;; http://blog.sciencenet.cn/u/lixujeremy
;; http://lixuworld.blogspot.com/

pro MVC_TinyMomory, fullnames, stroutpath
  ;; Create a blank Matrix 
  ok=Query_Tiff(fullnames[0], info)
  Cols=info.dimensions[0]
  Rows=info.dimensions[1]
  N_Dim=size(fullnames, /n_elements)
  OutMat=DblArr(Cols, Rows)  
  for ii=0, N_Dim-1 do Begin
    image=Read_Tiff(fullnames[ii], GEOTIFF=GEOVAR)  
    OutMat=OutMat>image  
  endfor
  
  ;;Write a result tiff
  Write_Tiff, stroutpath, OutMat, GEOTIFF=GEOVAR, /float
  print, stroutpath+' done!'


end

Matlab版

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
% Created by LI Xu
% Version 1.0
% June 21, 2016

% Description:
% MVC Matlab Version

% If you have any question about this code,
% please do not hesitate to contact me via E-mail: 
% jeremy456@163.com

% Blog:
% http://blog.sciencenet.cn/u/lixujeremy
% http://lixuworld.blogspot.com/

% Note;
% fullpathes: full file pathes in cell
% otpath:   output file path               

function MVC4Image(fullpathes, otpath)
    % First Image
    filepath=fullpathes{1};
    [output, geo]=geotiffread(filepath);
    info=geotiffinfo(filepath);
    [row, col]=size(output);
    
    output=output(:);
    for ii=2:length(fullpathes)
        filepath=fullpathes{ii};
        image=imread(filepath);
        image=image(:);        
        output=[output, image];
        output=max(output, [], 2);   
    end
        
    output=reshape(output, row, col);
    geotiffwrite(otpath, output, geo, 'GeoKeyDirectoryTag', info.GeoTIFFTags.GeoKeyDirectoryTag);
       
end

References

[1] Maximum-value composite procedure, from Wikipedia, the free encyclopedia.
Posted by at No comments:
Labels: ,

Thursday, March 10, 2016

ENVI+IDL+Matlab: Pixel Aggregate without Average

Introduction

在Fine Resolution向Coarse Resolution重采样过程中，ENVI已有的重采样算法Pixel Aggregate是将临近像元数值平均后得到输出像元数值。
Fig. 1
但，有时降采样的目的是要将Fine Resolution像元数值合并到Coarse Resolution像元数值，即后者是前者的累计之和，这种情况多见于人口密度分布等情况。
举例：Fig. 1左边Fine Resolution像元尺寸4×4 m，降采样至右边Coarse Resolution像元尺寸6×6 m，要求是合并输出新数值。如Fig. 2所示，以输出像元(2, 1)为例，说明处理过程，IDL表示像元在矩阵中的位置以列先行后为序，均起始于0，它由如下像元而得：
100%的输入(3, 2)=15；
50%的输入(3, 1)=9；
50%的输入(4, 2)=16；
25%的输入(4, 1)=10；
则该输出像元(2, 1)是由各部分合并之和得到的，如下：
1×15+0.5×9+0.5×16+0.25*10=30
Fig. 2
如上所述，全部输出像元数值：
(0, 0)：0×1+0.5×1+0.5×6+0.25×7=5.25
(1, 0)：0.5×1+1×2+0.5×8+0.25×7=8.25
(2, 0)：1×3+0.5×4+0.5×9+0.25×10=12
(3, 0)：1×5+0.5×4+0.5×11+0.25×10=15
(0, 1)：12+0.5×6+0.5×13+0.25×7=23.25
(1, 1)：14+0.5×8+0.5×13+0.25×7=26.25
(2, 1)：15+0.5×9+0.5×16+0.25×10=30
(3, 1)：17+0.5×16+0.5×11+0.25×10=33
(0, 2)：18+0.5×24+0.5×19+0.25×25=45.75
(1, 2)：20+0.5×19+0.5×26+0.25×25=48.75
(2, 2)：21+0.5×22+0.5×27+0.25×28=52.5
(3, 2)：23+0.5×22+0.5×29+0.25×28=55.5
(0, 3)：30+0.5×24+0.5×31+0.25×25=63.75
(1, 3)：32+0.5×31+0.5×26+0.25×25=66.75
(2, 3)：33+0.5×27+0.5×34+0.25×28=70.5
(3, 3)：35+0.5×34+0.5×29+0.25×28=73.5
得到降采样后的矩阵如下：
Fig. 3
简述处理过程，首先，ENVI+IDL处理输入数据以Pixel Aggregate（Average）方法重采样，得到6×6 m空间分辨率的输出结果，但此时的输出结果是ENVI Standard格式；第二步，转换为Geotiff格式文件；第三步，参阅文献[1]，将栅格数据乘以转换因子，祛除Average换为Summation，这一因子是目标分辨率（Coarse Resolution）与原始分辨率（Fine Resolution）之比值，最终得到像元数值表示累计之和的结果。注意：此处仅仅适用于影像空间分辨率X/Y相等的情况，对于X/Y不等的情况暂不适用。
示例影像降采样前后的空间覆盖范围完全一致。
Fig. 4
附上代码。

References

[1] How does the Pixel Aggregate method work when resizing data with ENVI?
Posted by at No comments:
Labels: , ,

Wednesday, March 9, 2016

Code: Archiver Code

[1] 人口压力密度分布.
[2] GDP密度分布.
[3] MODIS NDVI切换为AVHRR NDVI代码.
[4] 牲畜密度分布.
[5] 典型草原1989~2014年NDVI数据.
[6] 典型草原日照时数的处理代码。
[7] 二级土地利用图转一级土地利用图处理代码。
[8] 土地利用变化检测.
[9] 土地利用类型转换分析.
[10] 草地覆盖度空间分布.
[11] 典型草原1989~2014年生长季NDVI数据.

博士论文数据

[1] 指标数据.
[2] 检测Geotiff格式数据数值区间，单波段或多波段，加强版【2023-08】.
[3] 土地利用类型结合行政区划提取数值.
[4] 年度评价代码及归纳评价代码.
[5] 年度评价之全局评价靶心度贡献率代码（纵向评价靶心度贡献率代码）.
[6] 拆分多波段文件.
[7] 矩阵分级代码.
[8] 等级面积比重代码.（面向归纳评价）
[9] 提取生态安全指数（年度、归纳）.
[10] 归纳评价生态安全指数贡献率（全局、纵向）.
[11] 归纳权重识别（全局、纵向）.
[12] 重心识别（年度、归纳）.
[13] 指标分量的线性趋势及其显著性.
[14] 分组，返回索引位置，不同分组依次位于元胞之中.
[15] 归纳评价得分之子系统分量占比（全局、纵向）.
[16] 归纳评价之指标分量一元二次曲线变化率的变化率.
[17] 年度评价之全局评价公式修正.
[18] 指标评估得分.

25指标

[1] 删除部分靶心系数.
[2] 年度指标分量贡献提取（全局、纵向）.
[3] 栅格年度指标分量（全局、纵向）.
[4] 栅格多年指标分量（全局、纵向）.
[5] 多年指标分量贡献提取（全局、纵向）.
[6] 逐年准则层指标分量贡献（全局、纵向）.
[7] 多年准则层指标分量贡献（全局、纵向）.
[8] 重命名文件，达到重新排序的目的.
[9] 多年全局视角的权重.
[10] 添加指标代码.
[11] 趋势百分比.
[12] 逐年栅格线性回归参数a和b，显著性是1.
[13] 多年栅格一元二次回归参数a/b/c.
[14] 逐年趋势预测情景.
[15] 修剪情景数值分布在0~1之间.
[16] 矩阵分级（6等级）.
[17] 等级比重6等级（逐年、多年）.
[18] 一元二次情景计算ABC.
[19] 识别起止年份.
[20] 移动平均求算多年合成归纳数值.
[21] 移动平均三子系统指数分量（全局、纵向）.
[22] 移动平均指标层指数分量（全局、纵向）.
[23] 识别（移动平均）三子系统分地区指数分量.
[24] 识别（移动平均）指标层分地区指数分量.
[23] 移动平均时序演化.
[24] 逐年经济演化情景.
[25] 移动平均经济演化情景.
[26] 逐年情景渲染Package.
[27] 移动平均情景渲染Package.
[28] 逐年多情景生态安全数值分区提取.
[29] 移动平均多情景生态安全数值分区提取.
[30] 敏感度系数（逐年、移动平均）.
[31] 移动平均的指标评估得分.
[32] 时序规律情景计算.
[33] 情景抽取态势指数.
[34] 情景抽取等级比重.
[35] 情景渲染和构图成片.
[36] 斜率截距渲染和构图成片.
Posted by at No comments:
Labels: , , , ,

Friday, November 20, 2015

IDL: 删除影像中的波段

问题：删除影像中的若干波段，并生成新影像文件？
方法：IDL。练习数据中，test.tif包含13个波段，删除其中的1~N（小于13）个波段，得到结果result.tif。
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
;+
; :Author: LiXu
;-


;;Created by LI Xu
;;Version 1.0
;;19 October, 2014

pro Delete_Bands_Batch

  print, 'Time Begin:'+' '+systime()
  begintime=SYSTIME(1)
  
  ;;Source Directory
  SouDir='E:\Tools\rad\Rad'
  ;;Destination Directory
  DesDir='E:\Tools\rad\NEW'
  ;;The Bands to Delete
  band=[1]
  
  
  ;;Retrieve all files
  filespath=file_search(SouDir, '*.tif', count = num_files, /test_regular)
  for ii=0,  num_files-1 do begin
    imgpath=filespath(ii)
    ;;otImagepath
    filename=strsplit(imgpath, '\', count=str_count, /Extract)
    filename=filename(str_count-1)
    otpath=DesDir+'\'+filename
    
    Delete_Bands, imgpath, otpath, band
    
  endfor
  
  
  

  print, 'End!'
  print, 'Time End:'+' '+systime()
  endtime=systime(1)
  timespan=endtime-begintime
  print, 'Time Span:'+' '+string(timespan)+' s'







end


 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
;+
; :Author: LiXu
;-

;;Created by LI Xu
;;Version 1.0
;;18 October, 2014

;;Description:
;;Delete several bands from the input Image,
;;and generate a new Image.

pro Delete_Bands, infile, otfile, numoffile

  ;;Debug Block
  ;;infile='F:\Tools\Practice\test.tif'
  ;;otfile='F:\Tools\Practice\result.tif'
  ;;numoffile=[1]
  
  ;;Code Body
  imgfile=read_tiff(infile, GeoTIFF=GEOVAR)
  filesize=size(imgfile)
  Num_Bands=filesize(1)
  Bands=indgen(13)
  
  cols=filesize(2)
  rows=filesize(3)
  
  ;;The Bands to Delete
  Num_Dlet=N_Elements(numoffile)
  ;;The Bands to Save
  Num_Save=Num_Bands-Num_Dlet
  Remove, numoffile-1, Bands
  
  newfile=make_array(Num_Save, cols, rows, /float)
  for ii=0, Num_Save-1 do begin
    newfile[ii, *, *]=imgfile[Bands[ii], *, *]

  endfor

  ;;Write
  write_tiff, otfile, newfile, GEOTIFF=GEOVAR, /float





end


  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
pro remove,index, v1, v2, v3, v4, v5, v6, v7, v8, v9, v10, v11, v12, v13, v14, $
  v15, v16, v17, v18, v19, v20, v21, v22, v23, v24, v25
  ;+
  ; NAME:
  ;       REMOVE
  ; PURPOSE:
  ;       Contract a vector or up to 25 vectors by removing specified elements
  ; CALLING SEQUENCE:
  ;       REMOVE, index, v1,[ v2, v3, v4, v5, v6, ... v25]
  ; INPUTS:
  ;       INDEX - scalar or vector giving the index number of elements to
  ;               be removed from vectors.  Duplicate entries in index are
  ;               ignored.    An error will occur if one attempts to remove
  ;               all the elements of a vector.     REMOVE will return quietly
  ;               (no error message) if index is !NULL or undefined.
  ;
  ; INPUT-OUTPUT:
  ;       v1 - Vector or array.  Elements specifed by INDEX will be
  ;               removed from v1.  Upon return v1 will contain
  ;               N fewer elements, where N is the number of distinct values in
  ;               INDEX.
  ;
  ; OPTIONAL INPUT-OUTPUTS:
  ;       v2,v3,...v25 - additional vectors containing
  ;               the same number of elements as v1.  These will be
  ;               contracted in the same manner as v1.
  ;
  ; EXAMPLES:
  ;       (1) If INDEX = [2,4,6,4] and V = [1,3,4,3,2,5,7,3] then after the call
  ;
  ;               IDL> remove,index,v
  ;
  ;       V will contain the values [1,3,3,5,3]
  ;
  ;       (2) Suppose one has a wavelength vector W, and three associated flux
  ;       vectors F1, F2, and F3.    Remove all points where a quality vector,
  ;       EPS is negative
  ;
  ;               IDL> bad = where( EPS LT 0, Nbad)
  ;               IDL> if Nbad GT 0 then remove, bad, w, f1, f2, f3
  ;
  ; METHOD:
  ;       If more than one element is to be removed, then HISTOGRAM is used
  ;       to generate a 'keep' subscripting vector.    To minimize the length of
  ;       the subscripting vector, it is only computed between the minimum and
  ;       maximum values of the index.   Therefore, the slowest case of REMOVE
  ;       is when both the first and last element are removed.
  ;
  ; REVISION HISTORY:
  ;       Written W. Landsman        ST Systems Co.       April 28, 1988
  ;       Cleaned up code          W. Landsman            September, 1992
  ;       Major rewrite for improved speed   W. Landsman    April 2000
  ;       Accept up to 25 variables, use SCOPE_VARFETCH internally
  ;              W. Landsman   Feb 2010
  ;       Fix occasional integer overflow problem  V. Geers  Feb 2011
  ;       Quietly return if index is !null or undefined W.L. Aug 2011
  ;
  ;-
  On_error,2
  compile_opt idl2,strictarrsubs

  npar = N_params()
  nvar = npar-1
  if npar LT 2 then begin
    print,'Syntax - remove, index, v1, [v2, v3, v4,..., v25]'
    return
  endif

  if N_elements(index) EQ 0 then return

  vv = 'v' + strtrim( indgen(nvar)+1, 2)
  npts = N_elements(v1)

  max_index = max(index, MIN = min_index)

  if ( min_index LT 0 ) || (max_index GT npts-1) then message, $
    'ERROR - Index vector is out of range'

  if ( max_index Eq min_index ) then begin   ;Remove only 1 element?
    Ngood = 0
    if npts EQ 1 then message, $
      'ERROR - Cannot delete all elements from a vector'
  endif else begin


    ;  Begin case where more than 1 element is to be removed.   Use HISTOGRAM
    ;  to determine then indices to keep

    nhist = max_index - min_index +1

    hist = histogram( index)      ;Find unique index values to remove
    keep = where( hist EQ 0, Ngood ) + min_index

    if ngood EQ 0 then begin
      if ( npts LE nhist ) then message, $
        'ERROR - Cannot delete all elements from a vector'
    endif
  endelse

  imin = min_index - 1
  imax = max_index + 1
  i0 = (min_index EQ 0) + 2*(max_index EQ npts-1)
  case i0 of
    3: begin
      for i=0, nvar-1 do  $
        (SCOPE_VARFETCH(vv[i],LEVEL=0)) = $
        (SCOPE_VARFETCH(vv[i],LEVEL=0))[keep]
      return
    end

    1:  ii = Ngood EQ 0 ? imax + lindgen(npts-imax) : $
      [keep, imax + lindgen(npts-imax) ]
    2:  ii = Ngood EQ 0 ? lindgen(imin+1)               :  $
      [lindgen(imin+1), keep ]
    0:   ii = Ngood EQ 0 ? [lindgen(imin+1), imax + lindgen(npts-imax) ]  : $
      [lindgen(imin+1), keep, imax + lindgen(npts-imax) ]
  endcase

  for i=0,nvar-1 do  $
    (SCOPE_VARFETCH(vv[i],LEVEL=0)) =    $
    (SCOPE_VARFETCH(vv[i],LEVEL=0))[ii]

  return
end
Attachment.
Posted by at No comments:
Labels:

ENVI+IDL: 批量提取多幅影像共同区域的方法

问题：已知多幅裁剪之后的影像，他们存在着多一行（或列）的问题，将其中行（列）最小的影像作为多幅影像的共同区域，达到批量提取其他影像的目的？
方法：ENVI+IDL。
练习数据中，mask.tif是行（列）最小的影像，即为多幅影像的共同区域。1比mask.tif多一行（列），代码运行的结果是result，行（列）与mask.tif一致，其第一波段即是mask.tif，其他波段为1影像。
文件test.pro是代码的入口。代码运行结果与ENVI操作结果一致。
注意不同文件之间数据类型差别，可能在批处理过程中要适当修改代码中选择数据类型的规则（out_dt=max(in_dt)）。Referencing: IDL Data Types.
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
;;Created by LI Xu
;;Version 1.0
;;18 October, 2014

pro CommonRange_Batch

  print, 'Time Begin:'+' '+systime()
  begintime=SYSTIME(1)

  ;;Source Directory
  SouDir='E:\Tools\rad\Rad'
  ;;Destination Directory
  DesDir='E:\Tools\rad\NEW'
  ;;One Image
  imgpath='E:\Tools\rad\pro\mask.tif'
  
  
  ;;Retrieve all files
  filespath=file_search(SouDir, '*.tif', count=num_files, /test_regular)
  for ii=0, num_files-1 do begin
    file=filespath(ii)
    ;;file=strsplit(file, '.', /extract)
    ;;file=file(0)
    
    filename=strsplit(file, '\', /extract, count=strcount)
    filename=filename(strcount-1)
    otImage=DesDir+'\'+filename
    inImage=strarr(2)
    inImage[0]=imgpath
    inImage[1]=file
    
    CommonRange, inImage, otImage
    
    
  endfor
  

  print, 'End!'
  print, 'Time End:'+' '+systime()
  endtime=systime(1)
  timespan=endtime-begintime
  print, 'Time Span:'+' '+string(timespan)+' s'




end


 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
;+
; :Author: LiXu
;-

;;Created by LI Xu
;;Version 1.0
;;18 October, 2014

;;Description:
;;Make files in same Column and Row,
;;and upper left coodinate.
;;Reference:
;;http://www.exelisvis.com/Support/HelpArticlesDetail/TabId/219/ArtMID/900/ArticleID/4572/4572.aspx

pro CommonRange, imgpath, otimg

  ;;Debug Bolck
;  imgpath=strarr(2)
;  imgpath[0]='C:\AAAAA\practice\mask.tif'
;  imgpath[1]='C:\AAAAA\practice\1'
;  otimg='C:\AAAAA\practice\result'
  
  
  
  compile_opt idl2
  envi, /restore_base_save_files
  
  ;;Number of files
  num_files=N_Elements(imgpath)
  ;;open and gather file information in arrays
  in_fid=lonarr(num_files)
  in_nb=lonarr(num_files)
  in_dims=lonarr(5, num_files)
  in_dt=lonarr(num_files)
  
  for i=0L, num_files-1 do begin
    envi_open_file, imgpath[i], r_fid=r_fid
    if (r_fid eq -1) then begin
      print, imgpath[i]+' Open failed!'
      return
    endif

    envi_file_query, r_fid, ns=ns, nl=nl, nb=nb, dims=dims, data_type=dt
    in_fid[i]=r_fid
    in_nb[i]=nb
    in_dims[*,i]=dims
    in_dt[i]=dt
  endfor
  
  
  ;set up output fid, pos, and dims arrays
  out_fid = replicate(in_fid[0], in_nb[0])
  for i=1, num_files-1 do out_fid = [out_fid, replicate(in_fid[i], in_nb[i])]

  out_pos = lindgen(in_nb[0])
  for i = 1, num_files-1 do out_pos = [out_pos, lindgen(in_nb[i])]

  rep_dims = (intarr(in_nb[0])+1) # in_dims[*, 0]
  for i = 1, num_files-1 do $
    rep_dims = [rep_dims, (intarr(in_nb[i]) + 1) # in_dims[*, i]]
  out_dims = transpose(rep_dims)
  
  ;set the output projection and pixel size from the first file.
  ;save the result to disk and use max data type
  out_proj = envi_get_projection(fid=in_fid[0], pixel_size=out_ps)
  out_dt = min(in_dt)
  out_name=otimg
  
  ;call the layer stacking routine.
  ;Use nearest neighbor for the interpolation method.
  envi_doit, 'envi_layer_stacking_doit', fid=out_fid, pos=out_pos, dims=out_dims, $
    out_dt=out_dt, out_name=out_name, interp=0, out_ps=out_ps, $
    out_proj=out_proj, r_fid=r_fid,  /EXCLUSIVE

 

  
  
  
  
end
Attachment.
Posted by at No comments:
Labels: ,

Thursday, November 19, 2015

ENVI+IDL: Resampling

问题：已知一个影像，对另一影像像元分辨率按照已知影像进行重采样？
方法：ENVI+IDL。
练习数据包括两张GEOTIFF文件，具有相同的空间参考，STANDARD.tif是已知影像，像元分辨率：1000*1000METERS；reprojected.tif是被重采样影像，像元分辨率：799.232051*799.232051METERS。
（原来，我尝试Python+GDAL进行重投影和重采样，但没发现Python+GDAL重采样的代码，所以就分别做处理了。）
reprojected.tif与resampled.tif前后转换对比，如图 1（似乎没什么不同？）。
图 1
注意：土地利用数据的重采样方法应选择最邻近法，它基于离散数据在插值过程中不会产生新的数值，适用于类似土地利用分类数据的插值操作（NEAREST—Performs a nearest neighbor assignment, is the fastest of the interpolation methods. It is used primarily for discrete data, such as a land-use classification, since it will not change the values of the cells. The maximum spatial error will be one-half the cell size.）。 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
;;Created by LI Xu
;;Version 1.0
;;16 october, 2014


pro Resample_Doit_Batch
  print, 'Time Begin:'+' '+systime()
  begintime=SYSTIME(1)

  ;;Standard Image
  stdimgpath='E:\Tools\rad\pro\mask.tif'
  ;;Source Directory
  SouDir='E:\Tools\rad\Rad'
  ;;Destination Directory
  DesDir='E:\Tools\rad\NEW'
  ;;Interpolation method
  ;;Interp
  ;;0: Nearest neighbor
  ;;1: Bilinear
  ;;2: Cubic convolution
  ;;3: Pixel aggregate
  intermed=3
  
  ;;Retrieve all files
  files=file_search(SouDir, '*.tif', count=file_counts, /test_regular)
  for ii=0, file_counts-1 do begin
    file=files(ii)
    filename=strsplit(file, '\', count=str_num ,/EXTRACT)
    filename=filename(str_num-1)
    otfile=DesDir+'\'+filename
    ;print, otfile
    Resample_Doit, stdimgpath, file, intermed, otfile
  endfor
  
  print, 'End!'
  print, 'Time End:'+' '+systime()
  endtime=systime(1)
  timespan=endtime-begintime
  print, 'Time Span:'+' '+string(timespan)+' s'

end


 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
;;Created by LI Xu
;;Version 1.0
;;16 October, 2014

;;Resampling x/y resolution size
;;Reference:
;;http://www.cnblogs.com/myyouthlife/archive/2012/07/01/2571785.html

;;Modified by LI Xu
;;Version 1.1
;;October 19, 2015
;;Add the variable for interpolating method


pro Resample_Doit, standardtif, resampltif, intemed, ottif
  
  ;;standardtif='D:\ForLIU\MOD_NDVI.tif'
  ;;resampltif='D:\ForLIU\LI_reprojected.tif'
  ;;ottif='D:\ForLIU\LI_resampled.tif'
  
  ;;*********************************
  compile_opt IDL2
  envi, /restore_base_save_files
  
  ;;Open the tifs
  ;;MODELPIXELSCALETAG DOUBLE[3] XYZ
  stdimg=read_tiff(standardtif, GEOTIFF=GEOVAR)
  tarXsize=GEOVAR.MODELPIXELSCALETAG(0)
  tarYsize=GEOVAR.MODELPIXELSCALETAG(1)
  envi_open_file, resampltif, r_fid=fid_resampled
  resampltif=read_tiff(resampltif, GEOTIFF=GEOVAR)
  nowXsize=GEOVAR.MODELPIXELSCALETAG(0)
  nowYsize=GEOVAR.MODELPIXELSCALETAG(1)
  envi_file_query, fid_resampled, dims=dims, nb=nb
  ;;DIMS[0]: A pointer to an open ROI; use only in cases where ROIs define the spatial subset. Otherwise, set to -1L.
  ;;DIMS[1]: The starting sample number. The first x pixel is 0.
  ;;DIMS[2]: The ending sample number
  ;;DIMS[3]: The starting line number. The first y pixel is 0.
  ;;DIMS[4]: The ending line number
  pos = lindgen(nb)
  
  Xsize=tarXsize/nowXsize
  Ysize=tarYsize/nowYsize
  envi_doit, 'resize_doit', $
    fid=fid_resampled, pos=pos, dims=dims, $
    interp=intemed, rfact=[Xsize, Ysize], $
    out_name=ottif, r_fid=r_fid
  ;;Interp  
  ;;0: Nearest neighbor
  ;;1: Bilinear
  ;;2: Cubic convolution
  ;;3: Pixel aggregate

end

References

[1] How does the Pixel Aggregate method work when resizing data with ENVI ?
Posted by at No comments:
Labels: ,

Monday, November 16, 2015

ENVI+IDL：经纬坐标提取影像的DN值

野外定点工作，常常需要和影像数据相匹配，就是以野外采样数据与影像的DN值两者之间进行比较、分析。那么，野外的样点坐标信息是已知的，剩下的就是依照样地地理坐标提取影像上的对应DN值。
有两种方法可以实现提取DN值的目的，ENVI手动获取和ENVI+IDL编程提取。下面，分别介绍这两种方法的步骤和特点：

ENVI手动方法

打开影像，Load Band一下，在主图像窗口双击打开CursorLocation/Value窗口，动态显示鼠标在图像上游艺的位置。接下来，还是在主图像窗口右键打开Pixel Locator，如图 1设置经度、纬度，注意图 1中的Sample表示列，Line表示行，均从1开始计算。
图 1
设置好后单击左下角Apply，查看之前打开的CursorLocation/Value窗口，Data显示0.152500即为图 1对应经纬度像元的DN值。
图 2
ENVI手动方法比较简单，容易上手。但是对于多组坐标和多幅图像，这个方法就太吃力了。

ENVI+IDL提取DN值

提取思路：将经纬度坐标转换为图像上对应的行列号，最后提取之。
代码运行的准确性，以方法1中的坐标检验代码提取的DN与ENVI是否相同？
如图 3，红框中的数值与方法1提取结果一致，其他坐标也是一样，证明这个方法还是有效的。
图 3
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
;;Created by LI Xu
;;15 April, 2014
;;Version 1.0

;;Extraction based on the known GeoInformation(Lat&Long)
;;Description: csv, a fullpath of an .csv
;;                      inImage, a fullpath of an image
;;                      otDN, returned DNs

;;Reference Website: http://www.gisall.com/wordpress/?p=8319
;;Modification Record
;;22 October, 2014
;;Version 1.1 Add output the No(Column, Row)=Value in the Console  Column and Row start from 0



pro ExtratDN, csv, inImage, otDN
  ;;Debug Block
  ;;csv='F:\Test\Practice\DNVI\h27v06.csv'
  ;;inImage='F:\Test\Practice\DNVI\MOD13A3.2008.h27v06.tif'
  csv='F:\Test\Practice\Exercise\GeoInfo.csv'
  inImage='F:\Test\Practice\Exercise\MOD_2011_EVI.tif'
  
  
  ;;Retrieve long&lat from .csv
  txt=Read_Csv(csv, header=header)
  no=txt.field1
  lon=txt.field2 ;;Longtitude
  lat=txt.field3  ;;Latitude 
  
  ;;Convert to the Projection as same as the Imput Image
  ;1. Retrieve the Pojection Info from Input Image
  ENVI_OPEN_FILE, inImage, r_fid=fid
  ProjInfo=ENVI_GET_PROJECTION(fid=fid)
  ;2. Original projection of Lon/Lat
  OriProj=ENVI_PROJ_CREATE(/geographic)
  ;3. Process of Convert
  ENVI_CONVERT_PROJECTION_COORDINATES, $
    lon, lat, OriProj, $
    otX, otY, ProjInfo
    
  ;;Extract the Col & Row for the input Geoinfo
  ;1. Geolocation of Upper Left Point (0, 0) of the Input Image
  mapinfo=ENVI_GET_MAP_INFO(fid=fid)
  ULlon=mapinfo.mc(2)
  ULlat=mapinfo.mc(3)
  ;2. X/Y Size of the Input Image
  Xsize=mapinfo.PS(0)
  Ysize=mapinfo.PS(1)
  ;3. Calculate the Number of Col & Row
  Coll=(otX-ULlon)/Xsize
  Row=(otY-ULlat)/Ysize
  Row=fix(abs(Row)) ;;Subscript directly from 0 
  Coll=fix(abs(coll))
  
  ;;Extraction of DN
  ;1. Number of DNs
  N_DN=N_Elements(Row)
  otDN=fltarr(N_DN)
  ;2. Retrieve the DN one by one
  tiff=Read_Tiff(inImage, GeoTIFF=GeoVar)
  for ii=0, N_DN-1 do begin
    otDN(ii)=tiff(Coll(ii), Row(ii)) 
    print, strtrim(string(no(ii)),2)+'('+ strtrim(string(Coll(ii)), 2)+','+strtrim(string(Row(ii)), 2)+'):='+strtrim(string(otDN(ii)), 2)
  endfor
  

  
end
Attachment.
Updated on May-16-2025.
Posted by at No comments:
Labels:
Subscribe to: Comments (Atom)