pro hai_asf_test_emccd_cdf

;; Initial setting
thm_init
loadct2,0

;; Directory of CDF files
dir='D:\asf\syo\2017\0722\'
fname='hai_asf_syo_201707220002_v01.cdf'

;; Load & get image data (pointer is used for speed-up!!)
cdf2tplot,file=dir+fname
get_data,'image_raw',ptr=p

;; Make & plot N-S & E-W keogram

dprint,'test to create keograms...'
store_data,'image_raw_keo_ns',data={x:*p.x,y:reform((*p.y)[*,127,*])},dlim={spec:1}
store_data,'image_raw_keo_ew',data={x:*p.x,y:reform((*p.y)[*,*,127])},dlim={spec:1}

tplot,'image_raw_keo_??'

dprint,'test for keogram done!!'
stop

;; create 2-D image in the image coordinate

dprint,'test to plot images...'

tt='2017-07-22/00:02:30'

ind=nn(*p.x,tt)
min=2000.
max=7000.

title='HAI SYO w/o Filter'

window,xs=600,ys=600
porg=!p
!p.background=0 & !p.color=255

for ii=ind[0],ind[0]+100 do begin

        xs=n_elements((*p.y)[ii,*,0])
        ys=n_elements((*p.y)[ii,0,*])

        plot,[0,0],[0,0],xrange=[0,xs-1],yrange=[0,ys-1],/iso,/xsty,/ysty,$
                        xtickformat='(A1)',ytickformat='(A1)',position=[0.125,0.075,0.875,0.825],$
                        color=0
        xyouts,0.5,0.95,title+'!C!C'+time_string((*p.x)[ii],precision=3),charsize=1.6,/normal,alignment=0.5

        tvimage,bytescale(reform((*p.y)[ii,*,*]),range=[min0,max0]),/keep_aspect,/over

wait,0.01
;stop

endfor

dprint,'test for plotting images done!!!'
stop

;; Create mapped images obtained by HAI at Syowa
;; Read mapping table
;; the format of the mapping tables is 256 x 256 x 15 (number of altitude)

dprint,'test for mapping images...'

cdfid=cdf_open(dir+fname)

cdf_varget,cdfid,'glat',glat_tmp
cdf_varget,cdfid,'glon',glon_tmp
cdf_varget,cdfid,'elevation_angle',elev

cdf_close,cdfid

;; only plot the data points with elevation angle above 5 deg.

ind_valid=where(abs(elev) gt 5.0)

alt=[80,85,90,95,100,105,110,115,120,125,130,150,200,250,300]

;; mapping the images

for ii=ind[0],ind[0]+4 do begin ; loop for images at different times

	for jj=0.,14 do begin ; loop for mapping tables at different altitudes

		dattmp=reform((*p.y)[ii,*,*])
		dat=dattmp[ind_valid]

		;; extract the mapping table for desired altitude.
		;; alt=[80,85,90,95,100,105,110,115,120,125,130,150,200,250,300]

		glat0=reform(glat_tmp[*,*,jj])
		glon0=reform(glon_tmp[*,*,jj])

		glat=reform(glat0[ind_valid])
		glon=reform(glon0[ind_valid])

		map_set,-69,40,scale=7e6,title=time_string((*p.x)[ii],precision=3)	+$
			' Alt: '+strcompress(string(alt[jj],format='(i3)'),/remove_all)+' km',charsize=1.4

		mapped=map_patch(dat,glon,glat,xstart=xst,ystart=yst,/triangulate)

		tv,bytescale(mapped,min=min,max=max),xst,yst

		map_grid,/label,charsize=1.4
		map_continents,/hires

	endfor

	stop

endfor

dprint,'test for mapping images done!!!'

!p=porg

end