爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 4243|回复: 2

[求助] 这个程序的思路是什么?

[复制链接]

新浪微博达人勋

发表于 2014-1-1 13:29:28 | 显示全部楼层 |阅读模式

登录后查看更多精彩内容~

您需要 登录 才可以下载或查看,没有帐号?立即注册 新浪微博登陆

x
这个程序是IDL编写的,我完全不懂这门语言。但其目的大致如下:对于一个降水场,设定一个临界值。该区域内大于该临界值的格点共同组成了一个降水主体,而在该降水主体内显然存在着不同的一个个不连续的降水个体。该程序的主要目的就是区分开这些降水个体。依次在mask矩阵中用1.2.3……标记,而非降水主体成员用0表示。希望各位大神帮帮忙,告诉小弟作者的编程思想,要是能够转成fortran,必感激不尽!!!!!!!!!代码如下::
;pro zero_in_on_maxes,grid1,grid2,thresh,critmass,ss1,ss2,ss3,ss4,entity_mask


pro isolate,grid,thresh,ie,area,iloc,jloc,maxval,mask

;IDL version of fortran subroutine ISOLAT (from McidAS NAWT command)
;Originally written by Andy Negri
;This may not be optimized at all for IDL!

;--- Searches array grid and defines entities by a threshold value.
;
;--- Input parameters are 'grid' (the array to be "contoured") and
;--- 'thresh' (the threshold value for an entity).
;
;--- Output is:
;---   ie (# of entities),
;---   area (# pixels in each entity),
;---   iloc,jloc (coord's of centroid)
;---   maxval (maximum value in each entity)
;---
;--- Additional output is the array 'mask', a mask of the original array
;--- with points less than thresh in grid set to zero in mask and
;--- points greater than or equal to thresh in grid set to the entity
;--- identifier (ie) in mask.

      s=size(grid)
      nx=s(1)
      ny=s(2)

      maxentities=800
      maxentitysize=8000        ;maximum size of an entity
      itogo=make_array(maxentitysize,/int,value=0)
      jtogo=itogo
      area=make_array(maxentities,/long,value=0)
      iloc=area & jloc=area & maxval=area
      isq=[-1,0,1]     ;square search area

      mask=make_array(nx,ny,/int,value=0)
      ie=0

;--- Search for entities - loop
      for jc=1,ny-2 do begin & for ic=1,nx-2 do begin   ;100 loop

      if mask(ic,jc) gt 0 or grid(ic,jc) lt thresh then goto, label100

;--- Brightness greater than threshold, increment cloud counter ie
      ie=ie+1    ;Found one!

      if ie ge maxentities then begin
       print,'Too many entities - do not save any more'
       ie=ie-1
       goto, done
       endif

;--- Starting index for fragment of cloud ie
      ntogo=0 ;was originally 1
      ndone=-1 ;was originally 0

;--- Place point's location in circular buffer
      itogo(0)=ic ;was originally itogo(1)
      jtogo(0)=jc ;same for this one

;--- Update histogram
      area(ie)=area(ie)+1
      iloc(ie)=iloc(ie)+ic
      jloc(ie)=jloc(ie)+jc
      maxval(ie)=grid(ic,jc)

;--- Update mask, marking point as found
      mask(ic,jc)=ie

addonpoints:

      if ndone ge ntogo then begin
       iloc(ie)=round(float(iloc(ie))/float(area(ie)))
       jloc(ie)=round(float(jloc(ie))/float(area(ie)))
       goto, label100
       endif

      ndone=ndone+1
      if ndone gt maxentitysize then ndone=0   ;was originally ndone=1

;--- Get center of next 8-point search
      i=itogo(ndone)
      j=jtogo(ndone)

;--- Search the 8 neighbors
      for l=0,2 do begin & for k=0,2 do begin     ;loop 50
      if k eq 1 and l eq 1 then goto, label50
      id=i+isq(k)
      jd=j+isq(l)
      if id lt 0 or id gt nx-1 then goto, label50
      if jd lt 0 or jd gt ny-1 then goto, label50
      if mask(id,jd) gt 0 or grid(id,jd) lt thresh then goto, label50

;--- Brightness greater than threshold so it's another point of cloud ie
       ntogo=ntogo+1
       if ntogo gt maxentitysize then ntogo=0    ;was originally ntogo=1

;--- Save location in circular buffer
       itogo(ntogo)=id
       jtogo(ntogo)=jd

;--- Update # of pixels
       area(ie)=area(ie)+1
       iloc(ie)=iloc(ie)+id
       jloc(ie)=jloc(ie)+jd
       maxval(ie)=max([maxval(ie),grid(id,jd)])

;--- Update mask, marking point as found
      mask(id,jd)= ie

label50:
      endfor & endfor

;--- Check that the area is large enough (> 2 pixels)
;--- If not, scratch out this entity
      if area(ie) lt 3 then begin
       area(ie)=0
       iloc(ie)=0
       jloc(ie)=0
       maxval(ie)=0
       for l=0,2 do begin & for k=0,2 do begin     ;loop 60
        id=i+isq(k)
        jd=j+isq(l)
        if id lt 0 or id gt nx-1 then goto, label60
        if jd lt 0 or jd gt ny-1 then goto, label60
        mask(id,jd)=0
label60:
        endfor & endfor
       ie=ie-1
       goto, label100
       endif

      goto, addonpoints

label100:
endfor & endfor

;Clean up - shrink array sizes down to where there is data
area=area(0:ie)
iloc=iloc(0:ie)
jloc=jloc(0:ie)
maxval=maxval(0:ie)

done:
end

pro getentity,mask,entID,ss

;Get the grid box including entity 'entID'

s=size(mask) & nx=s(1) & ny=s(2)
xx=indgen(nx)#replicate(1,ny)
yy=replicate(1,nx)#indgen(ny)
ssm=where(mask eq entID,nmask)
if nmask ge 1 then begin
  xstart=max([   0,min(xx(ssm))-1])
  ystart=max([   0,min(yy(ssm))-1])
  xend  =min([nx-1,max(xx(ssm))+1])
  yend  =min([ny-1,max(yy(ssm))+1])
  index=indgen(nx,ny)
  ss=index(xstart:xend,ystart:yend)
endif else ss=-1
end

pro zero_in_on_maxes,grid1,grid2,thresh,critmass,ss1,ss2,ss3,ss4,entity_mask

;Given two grids of data, this procedure zeros in on all of the local maxes
;and returns the subscripts, 'ss', of the boxes containing data surrounding
;the four largest local maxes.  Each box is defined by the area enclosing
;the isopleth given by 'thresh', where the value in the box can originate in
;either grid1 or grid2.

;Take the maximum value of grid1 and grid2
grid=grid1>grid2      

;Isolate entities with local maxima
;nentities = number of entities isolated
;area      = number of grid points in each entity
;iloc,jloc = indeces of centroid of each entity
;maxval    = maximum value inside each entity
;mask      = mask where each entity is marked with its identifying number,
;              all other gridpoints are 0
isolate,grid,thresh,nentities,area,iloc,jloc,maxval,mask

if nentities gt 0 then begin
  ;Get the integrated value inside each entity  (i.e., area * average value)
  intvalue=fltarr(nentities)
  for i=0,nentities-1 do intvalue(i)=total(grid(where(mask eq i+1)))
  ;Rank entities in order from large to small
  order=reverse(sort(intvalue))
  endif

;Get the boxes around the four largest entities
;Mimimum integrated size is 'critmass'
nentities=min([nentities,4])
ss=-1 & ss1=-1 & ss2=-1 & ss3=-1 & ss4=-1
n=0
for i=0,nentities-1 do begin
  entID=order(i)+1
  if intvalue(order(i)) ge critmass then begin
    getentity,mask,entID,ss
    if n eq 0 then ss1=ss
    if n eq 1 then ss2=ss
    if n eq 2 then ss3=ss
    if n eq 3 then ss4=ss
    n=n+1
    endif
endfor
nentities=n

;New code to write entity mask to file
entity_mask=mask & entity_mask[*,*]=0
for i=0,nentities-1 do begin
  entID=order(i)+1
  newID=i+1
  ss=where(mask eq entID)
  entity_mask(ss)=newID
  endfor
sz=size(entity_mask)

;Do a final check for repeated entities
if nentities eq 4 and ss4(0) eq ss3(0) then ss4=-1
if nentities ge 3 and ss3(0) eq ss2(0) then ss3=-1
if nentities ge 2 and ss2(0) eq ss1(0) then ss2=-1


done:
end


密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2014-1-8 19:31:04 来自手机 | 显示全部楼层
求回复!!!
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2016-3-21 21:03:06 | 显示全部楼层

回帖奖励 +2 金钱

原谅我混点钱,完全看不懂
密码修改失败请联系微信:mofangbao
您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表