IDL 二维数组/图像的Sen斜率实现

最新推荐文章于 2024-01-09 14:58:50 发布

终于搞副业啦

最新推荐文章于 2024-01-09 14:58:50 发布

阅读量679

点赞数

文章标签： idl

本文链接：https://blog.csdn.net/weixin_50079751/article/details/116598601

版权

参考文献《基于Sen+Mann-Kendall 的北京植被变化趋势分析》

针对二维数组/图像：

function Sen_slope,x1,y1,x2,y2
  s = (y2 - y1)/(x2 - x1)
  return,s
  
pro Sen_slope_image
  COMPILE_OPT idl2
  e = ENVI(/headless)
  tic
  ;--------------------------------------------------------------------------------
  ;打开图像，将所有图像存成band，合成一个raster
  ;打开一个band，记为basraster,获取spatialref
  ;raster命名为bandarray
  ;波段数为count
  ;uri为outpath1
  ;--------------------------------------------------------------------------------
  ;计算出的斜率uri为outpath2
  ;--------------------------------------------------------------------------------
  ;输出的斜率中值uri为outpath3
  ;--------------------------------------------------------------------------------
  
  ;计算斜率
  newcount = (count * (count-1))/2 ;计算新波段数
  slope = objarr(newcount) ;存放所有计算出的斜率波段

  bandnum = size(bandarray,/n_elements);计算合成新波段数组的大小
  raster = e.openraster(outpath1)
  m = 0 
  for j = 0, bandnum -1 do begin
    for i = 0, bandnum -1 do begin
      if i lt j then begin
        band1 = raster.getdata(bands = i)
        band2 = raster.getdata(bands = j)
        s = Sen_slope(i,band1,j,band2) ;计算斜率
        newraster = enviraster(s,uri = e.GetTemporaryFilename(), SPATIALREF = basraster.spatialref)
        newraster.save ;存成raster才能用envitask
        slope[m] = newraster
        m = m + 1
      endif
    endfor
  endfor
  Task = ENVITask("BuildBandStack")
  Task.SPATIAL_REFERENCE  = basraster.spatialRef
  Task.INPUT_RASTER = slope
  Task.OUTPUT_RASTER_URI = outpath2
  Task.Execute
  ;--------------------------------------------------------------------------------
  ;计算中值
  newraster = e.openraster(outpath2)
  slopeband  = newraster.getdata();打开所有波段
  medianslope = fltarr(columns,rows)

  for t = 0, columns-1 do begin ;列
    for n = 0, rows-1 do begin ;行
      medianslope[t,n] = median(slopeband[t,n,*]);列，行，波段数
    endfor
  endfor
  newfile = enviraster(medianslope,uri = outpath3,SPATIALREF = basraster.spatialref)
  newfile.save

  toc
end