Skip to content
Snippets Groups Projects
Commit ad37dc24 authored by David Hoese's avatar David Hoese
Browse files

added windpro lidar initial idl from nadia

parent 2bce1813
No related branches found
No related tags found
No related merge requests found
function compute_winds,s,status=status,snr_threshold=snr_threshold
status = 0b
t = mean(s.t)
t_ave = max(s.t)-min(s.t)
dims = size(s.ur,/dimensions)
if ((n_elements(dims) ne 2) and (dims[0] lt 3)) then return,-1
el = !pi*s.el/180.0
az = !pi*s.az/180.0
z = s.r*median(sin(el))
xhat = sin(az)*cos(el)
yhat = cos(az)*cos(el)
zhat = sin(el)
ur = s.ur
snr = s.snr
if (n_elements(snr_threshold) ne 1) then snr_threshold = 0.0
mean_snr = total(s.snr,1,/nan)/total(finite(s.snr),1)
u = replicate(!values.f_nan,dims[1])
v = replicate(!values.f_nan,dims[1])
w = replicate(!values.f_nan,dims[1])
residual = replicate(!values.f_nan,dims[1])
corr = replicate(!values.f_nan,dims[1])
; print,' Computing wind profile'
for i=0,dims[1]-1 do begin
ur1 = ur[*,i]
snr1 = snr[*,i]
index = where(snr1 ge snr_threshold,count)
if (count ge 4) then begin
ur1 = ur1[index]
xhat1 = xhat[index]
yhat1 = yhat[index]
zhat1 = zhat[index]
a = fltarr(3,3)
b = fltarr(3)
a[0,0] = total(xhat1^2,/nan)
a[1,0] = total(xhat1*yhat1,/nan)
a[2,0] = total(xhat1*zhat1,/nan)
a[0,1] = a[1,0]
a[1,1] = total(yhat1^2,/nan)
a[2,1] = total(yhat1*zhat1,/nan)
a[0,2] = a[2,0]
a[1,2] = a[2,1]
a[2,2] = total(zhat1^2,/nan)
b[0] = total(ur1*xhat1,/nan)
b[1] = total(ur1*yhat1,/nan)
b[2] = total(ur1*zhat1,/nan)
ainv = invert(a,/double)
condition = norm(a)*norm(ainv) ; Condition Number ?
if(condition lt 1.0e4) then begin
c = ainv#b
u[i] = c[0]
v[i] = c[1]
w[i] = c[2]
ur_fit = xhat1*u[i]+yhat1*v[i]+zhat1*w[i]
residual[i] = sqrt(mean((ur_fit-ur1)^2))
corr[i] = correlate(ur_fit,ur1)
endif
endif
endfor
; Compute windspeed and direction
wspd = sqrt(u^2+v^2)
wdir = 180.0*atan(u,v)/!pi + 180.0
status = 1b
return,{t:t,dt:t_ave,z:z,az:s.az,el:median(s.el),naz:uint(dims[0]),ur:ur,$
mean_snr:mean_snr,snr:snr,u:u,v:v,w:w,$
wspd:wspd,wdir:wdir,$
residual:residual,corr:corr}
end
pro DisplayColorBar,range,title=title,vertical=vertical,margin=margin,$
c_range=c_range,tickname=tickname,ticks=ticks,$
charsize=charsize,ticklen=ticklen,font=font
if (n_elements(c_range) ne 2) then c_range = [1,254]
ncolors = c_range[1]-c_range[0]
colors = c_range[0] + indgen(ncolors)
if (keyword_set(vertical)) then begin
boundary = [0.0,range[0],1.0,range[1]]
dx = 1.0
dy = (range[1]-range[0])/ncolors
x = replicate(0.5,ncolors)
y = range[0] + dy*findgen(ncolors) + dy/2.0
if (keyword_set(margin)) then begin
xmargin = margin
ymargin = [0,0]
endif
endif else begin
boundary = [range[0],0.0,range[1],1.0]
dx = (range[1]-range[0])/ncolors
dy = 1.0
x = range[0] + dx*findgen(ncolors) + dx/2.0
y = replicate(0.5,ncolors)
if (keyword_set(margin)) then begin
xmargin = [0,0]
ymargin = margin
endif
endelse
if (n_elements(charsize) eq 0) then charsize=1.0
plot,[0],[0],xstyle=5,ystyle=5,xrange=boundary[[0,2]],$
yrange=boundary[[1,3]],xmargin=xmargin,$
ymargin=ymargin,charsize=charsize,/nodata,/noerase
xvertices = dx*[-1.0,1.0,1.0,-1.0,-1.0]/2.0
yvertices = dy*[-1.0,-1.0,1.0,1.0,-1.0]/2.0
for i=0,ncolors-1 do $
polyfill,x[i]+xvertices,y[i]+yvertices,color=colors[i],$
clip=boundary,noclip=0,/data
if (keyword_set(vertical)) then begin
axis,1.0,0,YAXIS=1,ytitle=title,ystyle=1,font=font,$
xmargin=xmargin,yrange=range,ytickname=tickname,$
charsize=charsize,yticks=ticks,yticklen=ticklen,/DATA
endif else begin
axis,0,1.0,XAXIS=1,xtitle=title,xstyle=1,font=font,$
ymargin=ymargin,xrange=range,xtickname=tickname,$
charsize=charsize,xticks=ticks,xticklen=ticklen,/DATA
endelse
end
pro ImagePlot_old,x,y,u,boundary,xtitle=xtitle,ytitle=ytitle,scale_range=scale_range,$
charsize=charsize,xax=xax,yax=yax,noerase=noerase,dx=dx,dy=dy, $
missing=missing,range=range,points=points,xticks=xticks,yticks=yticks,$
xticklen=xticklen,yticklen=yticklen,font=font,xminor=xminor,nodata=nodata,$
charthick=charthick
if (n_elements(scale_range) ne 2) then scale_range = [1,254]
if (n_elements(range) ne 2) then begin
umin = min(u,max=umax)
range = [umin,umax]
endif
if (n_elements(dx) eq 0) then dx = median((x-shift(x,1))[1:*])
if (n_elements(dy) eq 0) then dy = median((y-shift(y,1))[1:*])
info = size(u)
x1 = x
y1 = y
u1 = u
if (info[0] eq 2) then begin
x1 = rebin(x,info[[1,2]])
y1 = transpose(rebin(y,info[[2,1]]))
endif
index = where(finite(u1),count)
if (count ne 0) then begin
x1 = x1[index]
y1 = y1[index]
u1 = u1[index]
endif
if (count eq 0) then nodata=1b
if (n_elements(missing) ne 0) then begin
index = where(u1 ne missing,count)
if (count ne 0) then begin
x1 = x1[index]
y1 = y1[index]
u1 = u1[index]
endif else nodata=1b
endif
n = n_elements(u1)
color = scale_range[0]+bytscl(u1,min=range[0],max=range[1],top = scale_range[1]-scale_range[0])
plot,[0.0],[0.0],xstyle=5,ystyle=5,xrange=boundary[[0,2]],yrange=boundary[[1,3]],$
/nodata,noerase=noerase
if (~keyword_set(nodata)) then begin
if (keyword_set(points)) then begin
plots,x1,y1,color=color,psym=3,clip=boundary,noclip=0
endif else begin
xpoly = dx*[-1.1,1.1,1.1,-1.1,-1.1]/2.0
ypoly = dy*[-1.1,-1.1,1.1,1.1,-1.1]/2.0
for i=0L,n-1 do begin
x2 = x1[i]+xpoly
y2 = y1[i]+ypoly
polyfill,x2,y2,color=color[i],clip=boundary,noclip=0
endfor
endelse
endif
plots,boundary[[0,2,2,0,0]],boundary[[1,1,3,3,1]]
IF (keyword_set(xax)) then axis,0,boundary[1],XAXIS=0,charsize=charsize,xminor=xminor,$
xrange=boundary[[0,2]],xtitle=xtitle,xstyle=1,$
xticks=xticks,xticklen=xticklen,font=font,charthick=charthick,/DATA
IF (keyword_set(yax)) then axis,boundary[0],0,YAXIS=0,charsize=charsize,$
yrange=boundary[[1,3]],ytitle=ytitle,ystyle=1,$
yticks=yticks,yticklen=yticklen,font=font,charthick=charthick,/DATA
end
pro ImagePlot,x,y,u,boundary,xtitle=xtitle,ytitle=ytitle,scale_range=scale_range,$
charsize=charsize,xax=xax,yax=yax,noerase=noerase,dx=dx,dy=dy, $
missing=missing,range=range,points=points,xticks=xticks,yticks=yticks,$
xticklen=xticklen,yticklen=yticklen,font=font,xminor=xminor,nodata=nodata,$
xtickunits=xtickunits,xtickformat=xtickformat,ytickformat=ytickformat,$
xtickinterval=xtickinterval,ytickinterval=ytickinterval,$
charthick=charthick
if (n_elements(scale_range) ne 2) then scale_range = [1,254]
if (n_elements(range) ne 2) then begin
umin = min(u,max=umax)
range = [umin,umax]
endif
if (n_elements(dx) eq 0) then dx = median((x-shift(x,1))[1:*])
if (n_elements(dy) eq 0) then dy = median((y-shift(y,1))[1:*])
info = size(u)
x1 = x
y1 = y
u1 = u
if (info[0] eq 2) then begin
x1 = rebin(x,info[[1,2]])
y1 = transpose(rebin(y,info[[2,1]]))
endif
index = where(finite(u1),count)
if (count ne 0) then begin
x1 = x1[index]
y1 = y1[index]
u1 = u1[index]
endif
count = 1
if (n_elements(missing) ne 0) then begin
index = where(u1 ne missing,count)
if (count ne 0) then begin
x1 = x1[index]
y1 = y1[index]
u1 = u1[index]
endif
endif
n = n_elements(u1)
; u1 = (u1 > range[0]) < range[1]
; color = scale_range[0] + (scale_range[1]-scale_range[0])*(u1-range[0])/(range[1]-range[0])
; color = byte(color)
color = scale_range[0]+bytscl(u1,min=range[0],max=range[1],top = scale_range[1]-scale_range[0])
plot,[0.0],[0.0],xstyle=5,ystyle=5,xrange=boundary[[0,2]],yrange=boundary[[1,3]],$
xtickunits=xtickunits,ytickunits=ytickunits,xtickformat=xtickformat,/nodata,noerase=noerase
if ((count ne 0) and ~keyword_set(nodata)) then begin
if (keyword_set(points)) then begin
plots,x1,y1,color=color,psym=3,clip=boundary,noclip=0
endif else begin
xpoly = dx*[-1.1,1.1,1.1,-1.1,-1.1]/2.0
ypoly = dy*[-1.1,-1.1,1.1,1.1,-1.1]/2.0
for i=0L,n-1 do begin
x2 = x1[i]+xpoly
y2 = y1[i]+ypoly
polyfill,x2,y2,color=color[i],clip=boundary,noclip=0
endfor
endelse
endif
plots,boundary[[0,2,2,0,0]],boundary[[1,1,3,3,1]]
IF (keyword_set(xax)) then axis,0,boundary[1],XAXIS=0,charsize=charsize,xminor=xminor,$
xrange=boundary[[0,2]],xtitle=xtitle,xstyle=1,$
xticks=xticks,xticklen=xticklen,font=font,charthick=charthick,$
xtickunits=xtickunits,xtickformat=xtickformat,$
xtickinterval=xtickinterval,ytickinterval=ytickinterval,/DATA
IF (keyword_set(yax)) then axis,boundary[0],0,YAXIS=0,charsize=charsize,$
yrange=boundary[[1,3]],ytitle=ytitle,ystyle=1,$
yticks=yticks,yticklen=yticklen,font=font,charthick=charthick,$
ytickunits=ytickunits,ytickformat=ytickformat,$
xtickinterval=xtickinterval,ytickinterval=ytickinterval,/DATA
end
pro plot_vectors,x,y,wdir,boundary,scale,stride,missing
dims = size(wdir,/dimensions)
dir = wdir
x2d = rebin(x,dims)
y2d = transpose(rebin(y,dims[[1,0]]))
i = stride[0]*indgen(dims[0]/stride[0])
x2d = x2d[i,*]
y2d = y2d[i,*]
dir = dir[i,*]
j = stride[1]*indgen(dims[1]/stride[1])
x2d = x2d[*,j]
y2d = y2d[*,j]
dir = dir[*,j]
index = where((dir ne missing) and (x2d ge boundary[0]) and (x2d le boundary[2]) and $
(y2d ge boundary[1]) and (y2d le boundary[3]),npts)
if (npts eq 0) then return
x2d = x2d[index]
y2d = y2d[index]
dir = dir[index]*!pi/180.0 - !pi
for i=0L,npts-1 do begin
D = CONVERT_COORD(x2d[i],y2d[i],/TO_DEVICE)
ARROW,d[0],d[1],d[0]+scale*sin(dir[i]),d[1]+scale*cos(dir[i]),hsize=-0.3
endfor
end
pro plot_wspd_wdir_time_height,t1,z1,wspd,wdir,datestr,snr_threshold,wspd_scale,stride,vector_scale,zmax,png_path
basename = 'Halo_winds.'+datestr
year = long(strmid(datestr,0,4))
month = long(strmid(datestr,4,2))
day = long(strmid(datestr,6,2))
t = 24.0*(julday(1,1,1970,0,0,t1)-julday(month,day,year,0))
if (n_elements(t) le 1) then return
dt = median((t-shift(t,1))[1:*])
title = 'Halo winds during XPIA ' + datestr
snr_str = 'SNR > '+strtrim(string(snr_threshold,format='(F12.3)'),2)
boundary = [0.0,0.0,24.0,zmax/1000.0]
z = z1/1000.0
dims = size(wspd,/dimensions)
set_plot,'z'
SetColorTable,41
device,set_resolution=[700,500],set_character_size=[8,12],$
set_colors=256L
!p.color = 0
!x.margin = [8,4]
!y.margin = [0,0]
erase,255
!p.region = [0.9,0.2,1.0,0.8]
DisplayColorBar,wspd_scale,title='Wind Speed (ms!U-1!N)',margin=[0,6],charsize=1.2,ticklen=-0.2,/vertical
!p.region = [0.0,0.1,0.93,0.94]
imageplot,t,z,wspd,boundary,ytitle='Height (km)',scale_range=[1,254],charsize=1.2,$
/xax,/yax,range=wspd_scale,xtitle='Hour (UTC)',yticklen=-0.01,/noerase
plot_vectors,t,z,wdir,boundary,vector_scale,stride,-99999.0
xyouts,0.12,0.88,snr_str,charsize=1.2,/normal,alignment=0.0
xyouts,0.5,0.96,title,charsize=1.5,/normal,alignment=0.5
xyouts,0.99,0.01,'RK Newsom, '+systime(0),/normal,alignment=1.0,charsize=0.7
if (~file_test(png_path,/directory)) then file_mkdir,png_path
outfile = png_path+path_sep()+basename+'.png'
saveimage,outfile
end
function compute_winds,s,status=status,snr_threshold=snr_threshold
status = 0b
t = mean(s.t)
t_ave = max(s.t)-min(s.t)
dims = size(s.ur,/dimensions)
if ((n_elements(dims) ne 2) and (dims[0] lt 3)) then return,-1
el = !pi*s.el/180.0
az = !pi*s.az/180.0
z = s.r*median(sin(el))
xhat = sin(az)*cos(el)
yhat = cos(az)*cos(el)
zhat = sin(el)
ur = s.ur
snr = s.snr
if (n_elements(snr_threshold) ne 1) then snr_threshold = 0.0
mean_snr = total(s.snr,1,/nan)/total(finite(s.snr),1)
u = replicate(!values.f_nan,dims[1])
v = replicate(!values.f_nan,dims[1])
w = replicate(!values.f_nan,dims[1])
residual = replicate(!values.f_nan,dims[1])
corr = replicate(!values.f_nan,dims[1])
; print,' Computing wind profile'
for i=0,dims[1]-1 do begin
ur1 = ur[*,i]
snr1 = snr[*,i]
index = where(snr1 ge snr_threshold,count)
if (count ge 4) then begin
ur1 = ur1[index]
xhat1 = xhat[index]
yhat1 = yhat[index]
zhat1 = zhat[index]
a = fltarr(3,3)
b = fltarr(3)
a[0,0] = total(xhat1^2,/nan)
a[1,0] = total(xhat1*yhat1,/nan)
a[2,0] = total(xhat1*zhat1,/nan)
a[0,1] = a[1,0]
a[1,1] = total(yhat1^2,/nan)
a[2,1] = total(yhat1*zhat1,/nan)
a[0,2] = a[2,0]
a[1,2] = a[2,1]
a[2,2] = total(zhat1^2,/nan)
b[0] = total(ur1*xhat1,/nan)
b[1] = total(ur1*yhat1,/nan)
b[2] = total(ur1*zhat1,/nan)
ainv = invert(a,/double)
condition = norm(a)*norm(ainv) ; Condition Number ?
if(condition lt 1.0e4) then begin
c = ainv#b
u[i] = c[0]
v[i] = c[1]
w[i] = c[2]
ur_fit = xhat1*u[i]+yhat1*v[i]+zhat1*w[i]
residual[i] = sqrt(mean((ur_fit-ur1)^2))
corr[i] = correlate(ur_fit,ur1)
endif
endif
endfor
; Compute windspeed and direction
wspd = sqrt(u^2+v^2)
wdir = 180.0*atan(u,v)/!pi + 180.0
status = 1b
return,{t:t,dt:t_ave,z:z,az:s.az,el:median(s.el),naz:uint(dims[0]),ur:ur,$
mean_snr:mean_snr,snr:snr,u:u,v:v,w:w,$
wspd:wspd,wdir:wdir,$
residual:residual,corr:corr}
end
function get_ppi_data,file,ngates,home_point,min_gate=min_gate,status=status
status = 0b
s = read_hpl(file,min_gate=min_gate,/snr,status=status)
if (~status) then return,-1
status = 0b
ur = s.ur[*,0:ngates-1]
r = s.r[0:ngates-1]
az = (s.az + home_point) mod 360.0
snr = s.snr[*,0:ngates-1]
index = where(snr le 0.0,count)
if (count ne 0) then snr[index] = !values.f_nan
status = 1b
return,{t:s.t,r:r,ur:ur,snr:snr,az:az,el:s.el}
end
pro plot_vad_time_height,date,snr_threshold,stride,vector_scale,min_gate,wspd_scale,$
zmax,home_point,data_home,qlhome
yyyymm = strmid(date,0,6)
path = data_home + path_sep() + yyyymm + path_sep() + date
ngates = 100
files = file_search(path+path_sep()+'VAD_*.hpl',count=nfiles)
if (nfiles eq 0) then return
t = replicate(!values.f_nan,nfiles)
wspd = replicate(!values.f_nan,nfiles,ngates)
wdir = replicate(!values.f_nan,nfiles,ngates)
for i=0,nfiles-1 do begin
s = get_ppi_data(files[i],ngates,home_point,min_gate=min_gate,status=status)
if (status) then begin
vad = compute_winds(s,snr_threshold=snr_threshold,status=status)
t[i] = vad.t
z = vad.z
wspd[i,*] = vad.wspd
wdir[i,*] = vad.wdir
; help,vad,/struct
; stop
endif
endfor
index = where(finite(t),count)
if (count ne 0) then begin
t = t[index]
wspd = wspd[index,*]
wdir = wdir[index,*]
plot_wspd_wdir_time_height,t,z,wspd,wdir,date,snr_threshold,wspd_scale,stride,vector_scale,zmax,qlhome
endif
end
pro plot_wind_profiles,vad,outfile,snr_threshold,display_zmax,wspd_scale
caldat,julday(1,1,1970,0,0,vad.t),month,day,year,hour,minute,second
date_time = string(year,month,day,hour,minute,second,format='(I4.4,2I2.2,".",3I2.2)')
base = file_basename(outfile,'.png')
nbeams_str = strtrim(string(vad.naz),2)
snr_threshold_str = strtrim(string(snr_threshold,format='(F12.3)'),2)
el_str = strtrim(string(vad.el,format='(F10.1)'),2)+'!Uo!N'
dt_str = strtrim(string(vad.dt,format='(F10.1)'),2)+' s'
comment1 = base+', ' + date_time
comment2 = 'Nbeams = '+nbeams_str+', Elevation angle = ' + el_str+$
', scan duration = ' + dt_str+', SNR Threshold = ' + snr_threshold_str
z = vad.z
zmax = display_zmax
n = n_elements(z)
; Set up the graphics
set_plot,'z'
!p.color = 0
!p.background = 255
SetColorTable,41
!x.margin = [6,3]
!y.margin = [4,6]
device,set_colors=256l,set_character_size=[8,12],set_resolution=[800.0,400.0]
erase,255
A = FINDGEN(17) * (!PI*2/16.)
USERSYM, COS(A), SIN(A), /FILL
ytitle = '!3z (km AGL)'
!p.region = [0.0,0.0,0.4,1.0]
boundary = [0.0,0.0,360.0,zmax]
plot,[0],[0],xstyle=5,ystyle=1,xrange=boundary[[0,2]],yrange=boundary[[1,3]],$
ytitle=ytitle,charsize=1.0,/nodata,/noerase
plots,vad.wdir,z,color=254,thick=1,clip=boundary,noclip=0
plots,vad.wdir,z,psym=8,symsize=0.5,color=220,clip=boundary,noclip=0
axis,0.0,boundary[3],xaxis=1,xstyle=1,xrange=boundary[[0,2]],xtitle='Wind Direction (deg)',charsize=1.0,color=250
boundary = [wspd_scale[0],0.0,wspd_scale[1],zmax]
plot,[0],[0],xstyle=5,ystyle=1,xrange=boundary[[0,2]],yrange=boundary[[1,3]],$
ytitle=ytitle,charsize=1.0,/nodata,/noerase
plots,vad.wspd,z,color=0,thick=2,clip=boundary,noclip=0
plots,vad.wspd,z,psym=8,symsize=0.5,color=0,clip=boundary,noclip=0
axis,0.0,boundary[1],xaxis=0,xstyle=1,xrange=boundary[[0,2]],xtitle='Wind Speed (ms!U-1!N)',charsize=1.0
!p.region = [0.35,0.0,0.7,1.0]
boundary = [0.001,0.0,10.0,zmax]
plot,[0],[0],xstyle=5,ystyle=5,xrange=boundary[[0,2]],yrange=boundary[[1,3]],$
charsize=1.0,/nodata,/noerase,/xlog
plots,vad.mean_snr,z,color=0,thick=1,clip=boundary,noclip=0
plots,vad.mean_snr,z,color=0,thick=1,clip=boundary,noclip=0,psym=8,symsize=0.5
axis,1.0,boundary[1],xaxis=0,xstyle=1,xrange=boundary[[0,2]],xtitle='SNR',charsize=1.0,/xlog
boundary = [-10.0,0.0,10.0,zmax]
plot,[0],[0],xstyle=5,ystyle=5,xrange=boundary[[0,2]],yrange=boundary[[1,3]],charsize=1.0,/nodata,/noerase
plots,vad.w,z,color=254,thick=1,clip=boundary,noclip=0
plots,vad.w,z,color=220,thick=1,clip=boundary,noclip=0,psym=8,symsize=0.5
axis,0.0,boundary[3],xaxis=1,xstyle=1,xrange=boundary[[0,2]],xtitle='Vertical velocity (ms!U-1!N)',charsize=1.0,color=250
plots,[0.0,0.0],boundary[[1,3]],color=254,linestyle=1
plots,boundary[[0,2,2,0,0]],boundary[[1,1,3,3,1]]
!p.region = [0.65,0.0,1.0,1.0]
boundary = [0.0,0.0,10.0,zmax]
plot,[0],[0],xstyle=5,ystyle=5,xrange=boundary[[0,2]],yrange=boundary[[1,3]],charsize=1.0,/nodata,/noerase
plots,vad.residual,z,color=0,thick=1,clip=boundary,noclip=0
plots,vad.residual,z,color=0,thick=1,clip=boundary,noclip=0,psym=8,symsize=0.5
axis,0.0,boundary[1],xaxis=0,xstyle=1,xrange=boundary[[0,2]],xtitle='Residual (ms!U-1!N)',charsize=1.0
boundary = [0.0,0.0,1.0,zmax]
plot,[0],[0],xstyle=5,ystyle=5,xrange=boundary[[0,2]],yrange=boundary[[1,3]],charsize=1.0,/nodata,/noerase
plots,vad.corr,z,color=254,thick=1,clip=boundary,noclip=0
plots,vad.corr,z,color=220,thick=1,clip=boundary,noclip=0,psym=8,symsize=0.5
plots,boundary[[0,2,2,0,0]],boundary[[1,1,3,3,1]]
axis,0.0,boundary[3],xaxis=1,xstyle=1,xrange=boundary[[0,2]],xtitle='Correlation',charsize=1.0,color=250
xyouts,0.5,0.96,comment1,/normal,alignment=0.5,charsize=1.0
xyouts,0.5,0.92,comment2,/normal,alignment=0.5,charsize=1.0
xyouts,0.99,0.01,'Plot created on '+systime(0),/normal,alignment=1.0,charsize=0.7
path = file_dirname(outfile)
if (~file_test(path,/directory)) then file_mkdir,path
saveimage,outfile,/quiet
end
function get_ppi_data,file,min_gate=min_gate,rmax=rmax,zmax=zmax,status=status
status = 0b
parts = strsplit(file_basename(file),'_',/extract)
if (parts[0] ne 'VAD') then begin
print,'Invalid file'
return,-1
endif
s = read_hpl(file,/kilometers,status=status,rmax=rmax,/snr,min_gate=min_gate)
if (~status) then begin
print,'Error reading input file'
return,-1
endif
status = 0b
ur = s.ur
snr = s.snr
t = s.t
r = s.r
az = s.az
el = s.el
if (n_elements(zmax) ne 0) then begin
z = s.r*median(sin(s.el))
index = where(z le zmax,count)
if (count eq 0) then return,-1
r = r[index]
ur = ur[*,index]
snr = snr[*,index]
endif
status = 1b
return,{t:t,r:r,ur:ur,snr:snr,az:az,el:el}
end
pro write_data,vad,file,snr_threshold
wspd = vad.wspd
wdir = vad.wdir
w = vad.w
snr = vad.mean_snr
res = vad.residual
corr = vad.corr
missing = -9999.0
index = where(~finite(wspd),count)
if (count ne 0) then wspd[index] = missing
index = where(~finite(wdir),count)
if (count ne 0) then wdir[index] = missing
index = where(~finite(w),count)
if (count ne 0) then w[index] = missing
index = where(~finite(snr),count)
if (count ne 0) then snr[index] = missing
index = where(~finite(res),count)
if (count ne 0) then res[index] = missing
index = where(~finite(corr),count)
if (count ne 0) then corr[index] = missing
ngates = n_elements(vad.z)
caldat,julday(1,1,1970,0,0,vad.t),month,day,year,hour,minute,second
header = strarr(8)
header[0] = 'Date: '+string(year,month,day,format='(I4.4,2I2.2)')
header[1] = 'Time: '+string(hour,minute,second,format='(3I2.2)')
header[2] = 'SNR Threshold: '+strtrim(string(snr_threshold,format='(F10.3)'),2)
header[3] = 'Scan Duration: '+strtrim(string(vad.dt,format='(F10.2)'),2)+' s'
header[4] = 'Elevation angle: '+strtrim(string(vad.el,format='(F10.2)'),2)+' deg'
header[5] = 'Number of beams: '+strtrim(string(vad.naz,format='(I)'),2)+' deg'
header[6] = 'Number of gates: '+strtrim(string(ngates,format='(I)'),2)
header[7] = 'Missing Value: '+strtrim(string(missing,format='(F10.0)'),2)
path = file_dirname(file)
if (~file_test(path,/directory)) then file_mkdir,path
openw,unit,file,/get_lun
for i=0,n_elements(header)-1 do printf,unit,header[i]
printf,unit,'z(km) WindSpeed WindDir VertVel SNR Residual Correlation'
for i=0,ngates-1 do begin
c = strtrim(string(vad.z[i],wspd[i],wdir[i],w[i],snr[i],res[i],corr[i],$
format='(7F10.3)'),2)
printf,unit,c
endfor
free_lun,unit
end
pro process_vad_profile,files,snr_threshold,min_gate,rmax,zmax,wspd_scale,quicklook_path,out_data_path
nfiles = n_elements(files)
for i=0,nfiles-1 do begin
print,'Processing '+files[i]
s = get_ppi_data(files[i],min_gate=min_gate,rmax=rmax,zmax=zmax,status=status)
if (~status) then begin
print,'Error reading '+file
return
endif
vad = compute_winds(s,snr_threshold=snr_threshold,status=status)
if (~status) then begin
print,'Error computing VAD'
return
endif
; VAD_08_20150318_004018.hpl
base = file_basename(files[i],'.hpl')
date = (strsplit(base,'_',/extract))[2]
datafile = out_data_path + path_sep() + date + path_sep() + base + '.dat'
pngfile = quicklook_path + path_sep() + date + path_sep() + base + '.png'
plot_wind_profiles,vad,pngfile,snr_threshold,zmax,wspd_scale
write_data,vad,datafile,snr_threshold
endfor
end
function read_hpl,file,time_offset=time_offset,kilometers=kilometers,status=status,rmax=rmax,snr=snr,min_gate=min_gate
status = 0b
on_ioerror,bad
header = strarr(17)
openr,unit,file,/get_lun
readf,unit,header
gates = fix(strmid(header[2],16))
dr = float(strmid(header[3],22))
r = (findgen(gates) + 0.5)*dr
date_time = strtrim(strmid(header[9],11),2)
date = (strsplit(date_time,' ',/extract))[0]
date = strtrim(date,2)
year = 0 & month = 0 & day = 0
reads,date,year,month,day,format='(I4,2I2)'
data1 = fltarr(3)
data2 = fltarr(4,gates)
t = 0.0
az = 0.0
el = 0.0
readf,unit,t,az,el
readf,unit,data2
ur = reform(data2[1,*])
beta = reform(data2[3,*])
intensity = reform(data2[2,*])
while (~eof(unit)) do begin
readf,unit,data1
; print,data1
readf,unit,data2
t = [t,data1[0]]
az = [az,data1[1]]
el = [el,data1[2]]
ur = [[ur],[reform(data2[1,*])]]
beta = [[beta],[reform(data2[3,*])]]
intensity = [[intensity],[reform(data2[2,*])]]
endwhile
free_lun,unit
hour = fix(t)
index = where(hour eq 23,count)
if (count ne 0) then begin
indx = where(hour eq 0,cnt)
if (cnt ne 0) then t[indx] = t[indx] + 24.0
endif
if (n_elements(time_offset) eq 0) then time_offset = 0.0
jd = julday(month,day,year,t+time_offset)
t = 3600.0D*24.0D*(jd-julday(1,1,1970,0))
if (keyword_set(kilometers)) then r = r/1000.0
if (n_elements(rmax) eq 1) then begin
index = where(r le rmax,count)
if (count eq 0) then return,-1
r = r[index]
ur = ur[index,*]
beta = beta[index,*]
intensity = intensity[index,*]
endif
if (n_elements(min_gate) eq 1) then begin
if (min_gate ge n_elements(r)) then return,-1
r = r[min_gate:*]
ur = ur[min_gate:*,*]
beta = beta[min_gate:*,*]
intensity = intensity[min_gate:*,*]
endif
status = 1b
if (keyword_set(snr)) then begin
snr = intensity-1.0
return,{t:t,az:az,el:el,r:r,ur:transpose(ur),bscat:transpose(beta),snr:transpose(snr)}
endif
return,{t:t,az:az,el:el,r:r,ur:transpose(ur),bscat:transpose(beta),intensity:transpose(intensity)}
bad:
print,!err_string
close,/all
return,-1
end
data_home = 'C:\Users\rnewsom\projects\WindEnergy\XPIA\data\Halo\proc\dst_effected'
qlhome = 'ql\time_heights'
snr_threshold = 0.008
stride = [2,2]
vector_scale = 10.0
min_gate = 2
wspd_scale = [0.0,10.0]
zmax = 2000.0
home_point = 274.095
date = '20150405'
plot_vad_time_height,date,snr_threshold,stride,vector_scale,min_gate,wspd_scale,$
zmax,home_point,data_home,qlhome
end
\ No newline at end of file
quicklook_path = 'ql'
out_data_path = 'data'
snr_threshold = 0.008
min_gate = 2
rmax = 8.0
zmax = 2.0
wspd_scale = [0.0,20.0]
; specify a single file to process
;infile = 'C:\Users\rnewsom\projects\WindEnergy\XPIA\data\Halo\proc\dst_effected\201503\20150316\VAD_08_20150316_231851.hpl'
; or specify multiple files
data_home = 'C:\Users\rnewsom\projects\WindEnergy\XPIA\data\Halo\proc\dst_effected'
date = '20150316'
yyyymm = strmid(date,0,6)
path = data_home+path_sep()+yyyymm+path_sep()+date
infile = file_search(path+path_sep()+'VAD_08_'+date+'_*.hpl',count=nfiles)
process_vad_profile,infile,snr_threshold,min_gate,rmax,zmax,wspd_scale,quicklook_path,out_data_path
end
PRO SAVEIMAGE, FILE, BMP=BMP, PNG=PNG, PICT=PICT, JPEG=JPEG, TIFF=TIFF, GIF=GIF, $
QUALITY=QUALITY, DITHER=DITHER, CUBE=CUBE, QUIET=QUIET
;+
; NAME:
; SAVEIMAGE
;
; PURPOSE:
; Save the current graphics window to an output file (PNG by default).
;
; The output formats supported are:
; GIF 8-bit with color table,
; BMP 8-bit with color table,
; PNG 8-bit with color table,
; PICT 8-bit with color table,
; JPEG 24-bit true color,
; TIFF 24-bit true-color.
;
; Any conversions necessary to convert 8-bit or 24-bit images onscreen to
; 8-bit or 24-bit output files are done automatically.
;
; CATEGORY:
; Input/Output.
;
; CALLING SEQUENCE:
; SAVEIMAGE, FILE
;
; INPUTS:
; FILE Name of the output file (PNG format by default).
;
; OPTIONAL INPUTS:
; None.
;
; KEYWORD PARAMETERS:
; BMP Set this keyword to create BMP format (8-bit with color table).
; PNG Set this keyword to create PNG format (8-bit with color table).
; PICT Set this keyword to create PICT format (8-bit with color table).
; JPEG Set this keyword to create JPEG format (24-bit true color).
; TIFF Set this keyword to create TIFF format (24-bit true color).
; QUALITY If set to a named variable, specifies the quality for
; JPEG output (default 75). Ranges from 0 ("terrible") to
; 100 ("excellent"). Smaller quality values yield higher
; compression ratios and smaller output files.
; DITHER If set, dither the output image when creating 8-bit output
; which is read from a 24-bit display (default is no dithering).
; CUBE If set, use the color cube method to quantize colors when
; creating 8-bit output which is read from a 24-bit display
; (default is to use the statistical method). This may improve
; the accuracy of colors in the output image, especially white.
; QUIET Set this keyword to suppress the information message
; (default is to print an information message).
;
; OUTPUTS:
; None.
;
; OPTIONAL OUTPUTS:
; None
;
; COMMON BLOCKS:
; None
;
; SIDE EFFECTS:
; The output file is overwritten if it exists.
;
; RESTRICTIONS:
; Requires IDL 5.0 or higher (square bracket array syntax).
;
; EXAMPLE:
;
;openr, lun, filepath('hurric.dat', subdir='examples/data'), /get_lun
;image = bytarr(440, 330)
;readu, lun, image
;free_lun, lun
;loadct, 13
;tvscl, image
;saveimage, 'hurric.png'
;
; MODIFICATION HISTORY:
; Liam.Gumley@ssec.wisc.edu
; http://cimss.ssec.wisc.edu/~gumley
; $Id: saveimage.pro,v 1.1 2005/09/26 21:49:41 petty ds-idltools-idltools-1.28-0 $
;
; Copyright (C) 1999 Liam E. Gumley
;
; This program is free software; you can redistribute it and/or
; modify it under the terms of the GNU General Public License
; as published by the Free Software Foundation; either version 2
; of the License, or (at your option) any later version.
;
; This program is distributed in the hope that it will be useful,
; but WITHOUT ANY WARRANTY; without even the implied warranty of
; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
; GNU General Public License for more details.
;
; You should have received a copy of the GNU General Public License
; along with this program; if not, write to the Free Software
; Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
;-
rcs_id = '$Id: saveimage.pro,v 1.1 2005/09/26 21:49:41 petty ds-idltools-idltools-1.28-0 $'
;-------------------------------------------------------------------------------
;- CHECK INPUT
;-------------------------------------------------------------------------------
;- Check arguments
if (n_params() ne 1) then message, 'Usage: SAVEIMAGE, FILE'
if (n_elements(file) eq 0) then message, 'Argument FILE is undefined'
if (n_elements(file) gt 1) then message, 'Argument FILE must be a scalar string'
;- Check keywords
output = 'PNG'
if keyword_set(bmp) then output = 'BMP'
if keyword_set(gif) then output = 'GIF'
if keyword_set(pict) then output = 'PICT'
if keyword_set(jpeg) then output = 'JPEG'
if keyword_set(tiff) then output = 'TIFF'
if (n_elements(quality) eq 0) then quality = 75
;- Check for TVRD capable device
if ((!d.flags and 128)) eq 0 then message, 'Unsupported graphics device'
;- Check for open window
if (!d.flags and 256) ne 0 then begin
if (!d.window lt 0) then message, 'No graphics windows are open'
endif
;- Get display depth
depth = 8
if (!d.n_colors gt 256) then depth = 24
;-------------------------------------------------------------------------------
;- GET CONTENTS OF GRAPHICS WINDOW
;-------------------------------------------------------------------------------
version = float(!version.release)
;- Handle window devices (other than the Z buffer)
if (!d.flags and 256) ne 0 then begin
;- Copy the contents of the current display to a pixmap
current_window = !d.window
xsize = !d.x_size
ysize = !d.y_size
window, /free, /pixmap, xsize=xsize, ysize=ysize, retain=2
device, copy=[0, 0, xsize, ysize, 0, 0, current_window]
;- Set decomposed color mode for 24-bit displays
if (depth gt 8) then begin
if (version gt 5.1) then device, get_decomposed=entry_decomposed
device, decomposed=1
endif
endif
;- Read the pixmap contents into an array
if (depth gt 8) then begin
image = tvrd(order=0, true=1)
endif else begin
image = tvrd(order=0)
endelse
;- Handle window devices (other than the Z buffer)
if (!d.flags and 256) ne 0 then begin
;- Restore decomposed color mode for 24-bit displays
if (depth gt 8) then begin
if (version gt 5.1) then begin
device, decomposed=entry_decomposed
endif else begin
device, decomposed=0
if (keyword_set(quiet) eq 0) then $
print, 'Decomposed color was turned off'
endelse
endif
;- Delete the pixmap
wdelete, !d.window
wset, current_window
endif
;- Get the current color table
tvlct, r, g, b, /get
;- If an 8-bit image was read, reduce the number of colors
if (depth le 8) then begin
reduce_colors, image, index
r = r[index]
g = g[index]
b = b[index]
endif
;-------------------------------------------------------------------------------
;- WRITE OUTPUT FILE
;-------------------------------------------------------------------------------
case 1 of
;- Save the image in 8-bit output format
(output eq 'GIF') or (output eq 'BMP') or $
(output eq 'PICT') or (output eq 'PNG') : begin
if (depth gt 8) then begin
;- Convert 24-bit image to 8-bit
case keyword_set(cube) of
0 : image = color_quan(image, 1, r, g, b, colors=256, $
dither=keyword_set(dither))
1 : image = color_quan(image, 1, r, g, b, cube=6)
endcase
;- Sort the color table from darkest to brightest
table_sum = total([[long(r)], [long(g)], [long(b)]], 2)
table_index = sort(table_sum)
image_index = sort(table_index)
r = r[table_index]
g = g[table_index]
b = b[table_index]
oldimage = image
image[*] = image_index[temporary(oldimage)]
endif
;- In version 5.4, GIF output is not supported. Use PNG instead.
if(output eq 'GIF' and version ge 5.4) then begin
print,'GIF output not supported in this version -- using PNG instead'
output = 'PNG'
endif
;- Save the image
case output of
'GIF' : write_gif, file, image, r, g, b
'BMP' : write_bmp, file, image, r, g, b
'PNG' : write_png, file, image, r, g, b
'PICT' : write_pict, file, image, r, g, b
endcase
end
;- Save the image in 24-bit output format
(output eq 'JPEG') or (output eq 'TIFF') : begin
;- Convert 8-bit image to 24-bit
if (depth le 8) then begin
info = size(image)
nx = info[1]
ny = info[2]
true = bytarr(3, nx, ny)
true[0, *, *] = r[image]
true[1, *, *] = g[image]
true[2, *, *] = b[image]
image = temporary(true)
endif
;- If TIFF format output, reverse image top to bottom
if (output eq 'TIFF') then image = reverse(temporary(image), 3)
;- Write the image
case output of
'JPEG' : write_jpeg, file, image, true=1, quality=quality
'TIFF' : write_tiff, file, image, 1
endcase
end
endcase
;- Print information for the user
if (keyword_set(quiet) eq 0) then $
print, file, output, format='("Created ",a," in ",a," format")'
END
pro wab_color,red,green,blue
red = [ 0, 70,140,121,102, 82, 63, 47, 31, 15, $
0, 43, 86,130,173,216,255,255,255,255, $
255,216,201,186,170,155,140,163,186,209, $
232,245,255]
green =[255,255,255,216,178,140,102,140,178,216,$
255,255,255,255,255,255,255,235,216,197,$
178,140,132,124,117,109,102, 81, 61, 40,$
20, 10, 0]
blue = [ 0, 0, 0, 63,127,191,255,255,255,255,$
255,255,255,255,255,255,127, 95, 63, 31,$
0, 0, 12, 25, 38, 51, 63, 51, 38, 25,$
12, 6, 0]
x = (n_elements(red)-1)*findgen(!d.table_size)/(!d.table_size-1)
red = interpolate(red,x)
green = interpolate(green,x)
blue = interpolate(blue,x)
end
pro ddt_color,r,g,b
r = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,$
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,$
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,$
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,$
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,$
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,$
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, 15,$
29, 42, 56, 69, 82, 96,109,123,136,149,153,153,153,153,153,153,$
153,156,165,174,183,192,200,209,218,227,236,245,254,255,255,255,$
255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,$
255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,$
255,255,255,255,255,255,255,255,255,254,252,249,246,243,240,237,$
234,231,229,226,223,222,222,222,222,222,222,222,221,220,220,219,$
219,219,218,218,217,217,216,216,216,216,216,216,216,216,216,214,$
210,207,204,201,197,194,190,187,184,181,178,178,178,178,178,178,$
178,178,176,173,170,167,164,161,159,156,153,150,161,170,200]
g = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,$
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,$
0, 0, 0, 5, 16, 27, 38, 49, 60, 71, 82, 93,104,115,127,127,$
127,127,127,127,127,127,131,138,145,151,158,165,171,178,185,191,$
198,204,204,204,204,204,204,204,204,204,204,204,204,204,204,204,$
204,204,204,204,204,204,204,204,204,204,204,203,200,199,196,194,$
192,189,187,185,182,180,178,178,178,178,178,178,178,178,179,186,$
193,199,206,213,219,226,233,240,246,253,255,255,255,255,255,255,$
255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,$
255,255,255,255,250,243,237,230,223,217,211,204,198,191,184,180,$
180,180,180,180,180,180,180,172,162,153,144,134,125,116,107, 97,$
88, 79, 74, 74, 74, 74, 74, 74, 74, 73, 67, 60, 53, 47, 41, 34,$
27, 21, 15, 8, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,$
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,$
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,$
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 32, 0, 0]
b =[50,130,150,150,150,150,150,150,150,150,150,150,150,150,$
151,157,164,171,177,184,190,197,203,210,217,223,225,225,225,225,$
225,225,225,225,228,231,233,235,238,241,244,246,249,252,255,255,$
255,255,255,255,255,255,251,246,242,238,233,229,224,220,215,211,$
207,204,204,204,204,204,204,204,204,193,182,171,160,148,138,126,$
115,104, 93, 82, 76, 76, 76, 76, 76, 76, 76, 75, 68, 62, 55, 48,$
41, 35, 28, 21, 15, 8, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0,$
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,$
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,$
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,$
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,$
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,$
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, 22, 41, 59,$
78, 97,115,134,154,172,191,210,216,216,216,216,216,216,216,214,$
210,207,204,201,197,194,190,187,184,181,178,178,178,178,178,178,$
178,178,176,173,170,167,164,161,159,156,153,150,161,170,200]
x = (n_elements(r)-1)*findgen(!d.table_size)/(!d.table_size-1)
r = interpolate(r,x)
g = interpolate(g,x)
b = interpolate(b,x)
end
pro wind_direction,red,green,blue
; S W N E S
red = [255, 225, 0, 0, 255]
green = [ 0, 255, 0, 255, 0]
blue = [ 0, 0, 255, 0, 0]
x = (n_elements(red)-1)*findgen(!d.table_size)/(!d.table_size-1)
red = interpolate(red,x)
green = interpolate(green,x)
blue = interpolate(blue,x)
end
pro get_predefined_ct,n,r,g,b
loadct, (n > 0)
TVLCT,r,g,b,/GET
end
pro SetColorTable,CTNumber,ncolors=ncolors
case CTNumber of
41: wab_color,r,g,b
42: ddt_color,r,g,b
51: wind_direction,r,g,b
else: get_predefined_ct,CTNumber,r,g,b
endcase
n = !d.table_size-2
x = 1.0+n*findgen(n)/(n-1)
r = interpolate(r,x)
g = interpolate(g,x)
b = interpolate(b,x)
if (keyword_set(ncolors)) then begin
if ((ncolors gt 0) and (ncolors lt n)) then begin
xmax = max(x,min=xmin)
dx1 = (xmax-xmin)/ncolors
x1 = xmin + dx1*findgen(ncolors) + dx1/2.0
x2 = xmin + (xmax-xmin)*findgen(ncolors)/(ncolors-1)
r2 = interpolate(r,x2)
g2 = interpolate(g,x2)
b2 = interpolate(b,x2)
for i=0,ncolors-1 do begin
index = where(abs(x1[i]-x) le dx1/2.0,count)
if (count ne 0) then begin
r[index] = r2[i]
g[index] = g2[i]
b[index] = b2[i]
endif
endfor
endif
endif
r = [0,r,!d.table_size-1]
g = [0,g,!d.table_size-1]
b = [0,b,!d.table_size-1]
TVLCT,r,g,b
!P.COLOR = 0
!P.BACKGROUND = !d.table_size-1
end
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment