;===============================================================
;pattern_main_new4.pro
;Purpose: find the movement of smoke plume
;written by Jian Zeng
;2010-03-12
;===============================================================
pro find_move,target,search,glat,glon,outidx,$
target_xs=target_xs,target_ys=target_ys,mxcor=mxcor
pct=0.2
index=where(target gt 0.,nin)
outidx=[-1,-1]
mxcor=-999.0
if nin gt target_xs*target_ys*pct then begin
step=1
nx_sub=target_xs/step+1
ny_sub=target_ys/step+1
r=fltarr(nx_sub,ny_sub)-999.0
for m=0,nx_sub-1,step do begin
for n=0,ny_sub-1,step do begin
i1=m
i2=m+target_xs
j1=n
j2=n+target_ys
sub=search(i1:i2,j1:j2)
lat=glat(i1:i2,j1:j2)
lon=glon(i1:i2,j1:j2)
index2=where(lat gt -90 and lon gt -180. and sub gt 0.$
and target gt 0.,nin2)
if nin2 gt target_xs*target_ys*pct then begin
target_mean=mean(target(index2))
target_std=total((target(index2)-target_mean)^2)
sub_mean=mean(sub(index2))
sub_std=total((sub(index2)-sub_mean)^2)
up=total((sub(index2)-sub_mean)*(target(index2)-target_mean))
if sqrt(target_std*sub_std) gt 0. then $
r(m,n)=up/sqrt(target_std*sub_std)
endif
endfor
endfor
mxcor=max(r,location)
index=where(abs(mxcor-r) lt 1.e-2,cp)
if mxcor gt 0. then begin
outidx=array_indices(r,location)
endif
endif ;if nin gt target_xs*target_ys*pct
end
;----------------------------------------------------------------------------------------------------------
pro create_newimg,nx_g,ny_g,num,move,target_xs,target_ys,aod,newaod
newaod=fltarr(nx_g,ny_g)-9999.0
for inum=0,num-1 do begin
ii0=move(inum,0)
jj0=move(inum,1)
if ii0 gt 0 and jj0 gt 0 then begin
i1=(ii0-target_xs) > 0
i2=(ii0+target_xs) < (nx_g-1)
j1=(jj0-target_ys) > 0
j2=(jj0+target_ys) < (ny_g-1)
newaod(i1:i2,j1:j2)=aod(i1:i2,j1:j2)
endif
endfor
end
;------------------------------------------------------------------------------------
pro find_smokegrid,date_last,dir0,nx,ny,origx,origy,nx_g,ny_g,$
target_xs,target_ys,temp,nsmoke,data_suffix
infile_aod=file_search(dir0+'*'+date_last+data_suffix+'*',count=nfl_aod)
maxnum=nx*ny
temp=lonarr(maxnum,2)
nplumes=0
plumes=lonarr(maxnum,2)
nsmoke=0L
npix=target_xs*target_ys*0.1
if nfl_aod eq 1 then begin
saod=fltarr(nx_g,ny_g)
read_binary,infile_aod,saod
for inx=0,nx-1 do begin
for iny=0,ny-1 do begin
mini=origx(inx)
minj=origy(iny)
i1=(mini-target_xs/2) > 0
i2=(mini+target_xs/2) < (nx_g-1)
j1=(minj-target_ys/2) > 0
j2=(minj+target_ys/2) < (ny_g-1)
target=saod(i1:i2,j1:j2)
ind=where(target gt 0.,cp)
if cp gt npix then begin
temp(nsmoke,0)=mini
temp(nsmoke,1)=minj
nsmoke=nsmoke+1
endif
endfor
endfor
print,'Smoke Grids (#, max, min) at ',date_last,' = ',$
nsmoke,max(saod),min(saod(where(saod gt -9)))
endif
end
;=====================================================================================
;main program
;=====================================================================================
pro pattern_main_new4,file,instr,dir_smoke,yr,jday,time,yr_last,month_last,nday_last,$
dir_last,nx_g,ny_g,newaod,glat,glon,data_suffix,target_xs=target_xs,$
target_ys=target_ys,screen,maxang=maxang,minang=minang
if not keyword_set(target_xs) then target_xs=32
if not keyword_set(target_ys) then target_ys=32
nx=fix(nx_g/target_xs)-2L
ny=fix(ny_g/target_ys)-2L
if n_elements(file) lt 1 then begin
print, 'Pattern_main_new4.pro: No AOD file found'
stop
endif else begin
convert_julianday,yr,jday,mon,day
month=string(mon,format='(i2.2)')
nday=string(day,format='(i2.2)')
date0=yr+month+nday+time
print, 'starting from ',date0
pos=strpos(file(0),'_'+instr+'_US.all.aod.gz')
times=strmid(file,pos-13,11)
part=where(times eq yr+jday+time,nfl)
dir0=dir_smoke+yr+month+nday+'/'
if nfl ne 1 then print,'Pattern_main_new4.pro: no AOD file found at '+time+'Z' $
else if part(0) eq 0 then $
print,'no previous AOD file can be found' $
else begin
infile=file(part)
infile_last=file(part-1)
origx=indgen(nx)*target_xs+target_xs*3/2
origy=indgen(ny)*target_ys+target_ys*3/2
maxnum=nx*ny
move=lonarr(maxnum,2)-9999
corr=fltarr(maxnum)-9999.0
time_last=strmid(times(part-1),7,4)
date_last=yr+month+nday+time_last
if time_last gt time then begin
dir0=dir_last
date_last=yr_last+month_last+nday_last+time_last
endif
find_smokegrid,date_last,dir0,nx,ny,origx,origy,nx_g,ny_g,$
target_xs,target_ys,temp,nsmoke,data_suffix
if nsmoke eq 0L and part(0) gt 1 then begin
print, 'no smoke grid found 30 minutes ago, search one hour ago...'
;added on Feb. 11,2008
infile_last=file(part-2)
time_last=strmid(times(part-2),7,4)
if time_last gt time then begin
dir0=dir_last
date_last=yr_last+month_last+nday_last+time_last
endif else date_last=yr+month+nday+time_last
find_smokegrid,date_last,dir0,nx,ny,origx,origy,nx_g,ny_g,$
target_xs,target_ys,temp,nsmoke,data_suffix
endif ;if nsmoke eq 0L and part(0) ge 1
if nsmoke eq 0L then begin
print,'no smoke grid found within past 60 minutes'
endif else begin
print,'find pattern at -- ',date_last
aod0=read_GASP_sca_loose(infile,nx_g=nx_g,ny_g=ny_g,bsca=sca_ang0)
aod=read_GASP_sca_loose(infile_last,nx_g=nx_g,ny_g=ny_g,bsca=sca_ang)
if screen eq 1 then begin
if maxang ne '' then $
ind=where((sca_ang0 le minang and sca_ang0 gt -998.0) or sca_ang0 ge maxang,cp) $
else ind=where(sca_ang0 le minang and sca_ang0 gt -998.0,cp)
if cp gt 0 then aod0(ind)=-9999.0
if maxang ne '' then $
ind=where((sca_ang le minang and sca_ang gt -998.0) or sca_ang ge maxang,cp) $
else ind=where(sca_ang le minang and sca_ang gt -998.0,cp)
if cp gt 0 then aod(ind)=-9999.0
endif
ind=where(glat gt -90. and glon gt -180.)
for inum=0,nsmoke-1L do begin
mini=temp(inum,0)
minj=temp(inum,1)
if mini ge origx(0) and minj ge origy(0) and $
mini le origx(nx-1) and minj le origy(ny-1) then begin
i1=(mini-target_xs/2) > 0
i2=(mini+target_xs/2) < (nx_g-1)
j1=(minj-target_ys/2) > 0
j2=(minj+target_ys/2) < (ny_g-1)
target=aod0(i1:i2,j1:j2)
i1=(mini-target_xs) > 0
i2=(mini+target_xs) < (nx_g-1)
j1=(minj-target_ys) > 0
j2=(minj+target_ys) < (ny_g-1)
search=aod(i1:i2,j1:j2)
find_move,target,search,glat,glon,outidx,$
target_xs=target_xs,target_ys=target_ys,$
mxcor=mxcor
if outidx(0) ne -1 then begin
move(inum,0)=(outidx(0)-target_xs/2+mini) $
> (origx(0)-target_xs/2) < (origx(nx-1)+target_xs/2)
move(inum,1)=(outidx(1)-target_ys/2+minj)$
> (origy(0)-target_ys/2) < (origy(ny-1)+target_ys/2)
endif else begin
move(inum,0)=mini
move(inum,1)=minj
endelse
corr(inum)=mxcor
endif else begin ;if mini gt 0 and minj gt 0
move(inum,0)=mini
move(inum,1)=minj
endelse
endfor ;for inum=0,num-1L do begin
create_newimg,nx_g,ny_g,nsmoke,move,target_xs,target_ys,aod0,newaod
endelse ;if nsmoke eq 0L
endelse ;if nfl ne 1 or part(0) eq 0
endelse ;if nfl0 eq 0
end