; file: sdo_stx_fulldisk.pro ; init: Dec 28 2018 Rob Rutten Deil ; last: Jan 6 2019 Rob Rutten Deil ;+ pro sdo_stx_fulldisk,sdowav,stxname,$ rotate_up=rotate_up,markcenter=markcenter,marknorth=marknorth,$ xrange=xrange,yrange=yrange,thick=thick,verbose=verbose,$ alignfile=alignfile,$ fulldiskdir=fulldiskdir,level2dir=level2dir,resultdir=resultdir,$ datetime=datetime,solx=solx,soly=soly,nx_stx=nx_stx,ny_stx=ny_stx,$ px_stx=px_stx,angle_stx=angle_stx ; PURPOSE: write fits file with full-disk or partial SDO image ; with STX field of view outlined; optionally rotate to "up" ; ; CALL: ; see above ; ; MANDATORY INPUTS: ; sdowav: string specifying one of the following SDO "channels": ; AIA: '94', '131', '171', '193', '211', '304', '335', '1600', '1700' ; HMI: 'magnetogram', 'continuum' ; stxname: string specifying STX, eg 'sst' (same as in sdo_stx_align.pro) ; either alignfile (default sdo_stxname_align/sdo-stxname-align.dat) or if ; that is not specified or the default does not exist or if you ; want to override its values then specify all or a choice of: ; datetime: as '2013.12.21_12:00' or '2017-06-14T16:41:08.563" ; solx,soly: solar (X,Y) coordinates center STX field (arcsec) ; if datetime differs more than 5 min from the datetime in alignfile ; then you must re-specify solx, soly (use rot_xy.pro) ; nx_stx, ny_stx: STX image size in px_stx ; px_stx: STX pixel size in arcsec ; angle_stx: rotation angle NESW (CCW) of STX field from solar North ; ; OPTIONAL KEYWORD INPUTS: ; rotate_up = 1/0: rotate to nearest limb from STX ("up") (0) ; markcenter = 1/0: put black cross at image center (0) ; marknorth = 1/0: add black line from disk center to North (0) ; xrange=xrange: 2-elem X selection resulting image (full [0,4095]) ; yrange=yrange: 2-elem Y selection resulting image (full [0,4095]) ; thick=thick: thickness STX frame, center cross, North line (8) ; verbose = 1/0: more/less messaging ; fulldiskdir: main dir (default fulldisk) ; level2dir: store level "1.5" image (fulldiskdir/level2) ; resultdir: dir to store resulting image (fulldiskdir/annot) ; ; OUTPUTS: ; aia_prepped level2 image in level2dir for future calls same datetime ; imagefile = .fits file with full or partial image in resultdir ; ; HISTORY: ; Dec 28 2018 RR: start, mostly from sdo_featurelocator.pro ;- ; answer no-parameter query if (n_params(0) lt 2) then begin sp,sdo_stx_fulldisk return endif ; defaults for keywords if (n_elements(rotate_up) eq 0) then rotate_up=0 if (n_elements(markcenter) eq 0) then markcenter=0 if (n_elements(marknorth) eq 0) then marknorth=0 if (n_elements(xrange) eq 0) then xrange=[0,-1] if (n_elements(yrange) eq 0) then yrange=[0,-1] if (n_elements(thick) eq 0) then thick=8 if (n_elements(verbose) eq 0) then verbose=0 if (n_elements(alignfile) eq 0) then alignfile=$ 'sdo_'+stxname+'_align/sdo-'+stxname+'-align.dat' if (n_elements(fulldiskdir) eq 0) then fulldiskdir='fulldisk' if (n_elements(level2dir) eq 0) then level2dir=fulldiskdir+'/level2' if (n_elements(resultdir) eq 0) then resultdir=fulldiskdir+'/annot' if (n_elements(datetime) eq 0) then datetime='' if (n_elements(solx) eq 0) then solx=-1000 if (n_elements(soly) eq 0) then soly=-1000 if (n_elements(nx_stx) eq 0) then nx_stx=0 if (n_elements(ny_stx) eq 0) then ny_stx=0 if (n_elements(px_stx) eq 0) then px_stx=0 if (n_elements(angle_stx) eq 0) then angle_stx=-1000 ; define permitted SDO (AIA + HMI) request strings wavs=['94','131','171','193','211','304','335','1600','1700',$ ; AIA 'magnetogram','continuum'] ; HMI nwavs=n_elements(wavs) hmiseries=['hmi.M_45s','hmi.Ic_45s','hmi.V_45s'] ; check wav specification iwav=-1 for iw=0,nwavs-1 do if (wavs[iw] eq sdowav) then iwav=iw if (iwav eq -1 ) then begin print,' ##### sdo_stx_fulldisk abort: sdowav invalid = ',sdowav print,' ===== valid: ',wavs return endif ; make dirs if not yet existing spawn,'mkdir -p '+fulldiskdir spawn,'mkdir -p '+level2dir spawn,'mkdir -p '+resultdir ; =========== get location parameters from alignfile if it exists if (file_test(alignfile)) then begin ; read file values if (verbose) then print,' ----- sdo_stx_fulldisk reads '+alignfile openr,42,alignfile readf,42,angle_stx_file,px_stx_file,xasym,yasym readf,42,nt_stx,cadence_stx,nx_stx_file,ny_stx_file,nxrev,nyrev readf,42,solx_stx_file,soly_stx_file,angle_limb,muobs readf,42,itsample_stx,itsample_sdo datetime_file='' readf,42,datetime_file close,42 ; check for overrides in optional parameters if (datetime ne '') then begin diftime=abs(anytim2tai(datetime)-anytim2tai(datefime_file)) ; if file datetime within 5 min that is close enough to use file params if (diftime lt 300) then datetime=datetime_file else $ if (solx eq -1000 or soly eq -1000) then begin print,' ##### sdo_stx_fulldisk abort: '+$ 'datetime over 5 min from file value needs solx, soly' retall endif endif else datetime=datetime_file if (solx eq -1000) then solx=solx_stx_file if (soly eq -1000) then soly=soly_stx_file if (nx_stx eq 0) then nx_stx=nx_stx_file if (ny_stx eq 0) then ny_stx=ny_stx_file if (px_stx eq 0) then px_stx=px_stx_file if (angle_stx eq -1000) then angle_stx=angle_stx_file endif else begin if (nx_stx*ny_stx*px_stx eq 0 or solx eq -1000 or soly eq -1000 or $ angle_stx eq -1000 or datetime eq '') then begin print,' ##### sdo_stx_fulldisk abort: no align file and no parameters' retall endif endelse ; compute angle between North above and direction to nearest limb above ; (opposite to disk center; radial = projected vertical features up) angle_limb=limbangle(solx,soly) ; level2 file string for this datetime split=strmid(datetime,[0,5,8,11,14],[4,2,2,2,1]) ;RR f**k IDL datetimestr2=split[0]+split[1]+split[2]+'_'+split[3]+split[4] wavslevel2='*'+datetimestr2+'*'+$ ['_0094','_0131','_0171','_0193','_0211','_0304','_0335','_1600','_1700',$ '_mag','_cont']+'.fits' ; define output file wavsout=['aia94','aia131','aia171','aia193','aia211','aia304','aia335',$ 'aia1600','aia1700','hmimag','hmicont'] imagefile=wavsout[iwav]+'_'+datetimestr2+'_annot.fits' ; ================ get image ; get existing level2 image if present (within 10 minutes of datetime) if (file_test(level2dir+'/'+wavslevel2[iwav])) then begin infile=findfile(level2dir+'/'+wavslevel2[iwav]) ; take the first if multiple in this 10-min slot if (verbose) then print,' ----- sdo_stx_fulldisk uses existing '+infile[0] read_sdo,infile[0],indexim,im15,/uncomp_delete endif else begin if (verbose) then print,' ----- sdo_stx_fulldisk dowloads image '+wavs[iwav] t0=anytim2tai(datetime) if (iwav lt 9) then begin datetime1=anytim2utc(t0-20,/vms) datetime2=anytim2utc(t0+20,/vms) ssw_jsoc_time2data,datetime1,datetime2,ind1,im1,ds='aia.lev1',$ waves=[fix(wavs[iwav])],max_files=1,/uncomp_delete,/comp_delete end else begin datetime1=anytim2utc(t0+36,/vms) ; UTC lags TAI by Delta T approx 36 s datetime2=anytim2utc(t0+90,/vms) ssw_jsoc_time2data,datetime1,datetime2,ind1,im1,ds=hmiseries[iwav-9],$ max_files=1,/uncomp_delete,/comp_delete endelse ; aia_prep to "level 1.5" which I call level2 aia_prep,ind1,im1,indexim,im15,/do_write_fits,outdir=level2dir endelse ; ============= image manipulations ; rescale image intensity with default clips lev2wav=sdowav if (iwav eq 9) then lev2wav='mag' if (iwav eq 10) then lev2wav='cont' if (iwav eq 11) then lev2wav='dop' ; set RR standard clips (in sdo_intscale.pro) clipmin=-1 clipmax=-1 ; magnetogram: clip to see weaker fields if (iwav eq 9) then begin clipmin=-1000 clipmax=1000 endif ; intscale and histo_opt imscl=sdo_intscale(indexim,im15,lev2wav,clipmin=clipmin,clipmax=clipmax) if (iwav lt 9) then imscl=histo_opt_rr(imscl) ; bytscale or bytzero (magnetogram, Dopplergram)) if (iwav eq 9 or iwav eq 11) then $ image=bytzero(imscl) else image=bytscl(imscl) ; plot in invisible window window,22,/pixmap,xsize=4096,ysize=4096 tv,image ; optionally mark disk center with a black plus if (markcenter) then $ plots,2048.5,2048.5,psym=1,symsize=20,thick=thick,color=0,/device ; optionally add black line = half meridian from center to North if (marknorth) then $ plots,[2048.5,2048.5],[2048.5,4095],thick=thick,color=0,/device ; insert overlay frame showing tilted STX field of view nx=nx_stx*px_stx/0.6 ny=ny_stx*px_stx/0.6 xloc=fix(solx/0.6+2048.5) ; indexim.crpix = 2048.5 = some subtility? yloc=fix(soly/0.6+2048.5) llxy=rotatevec([-nx/2,-ny/2],-angle_stx)+[xloc,yloc] ulxy=rotatevec([-nx/2,+ny/2],-angle_stx)+[xloc,yloc] lrxy=rotatevec([+nx/2,-ny/2],-angle_stx)+[xloc,yloc] urxy=rotatevec([+nx/2,+ny/2],-angle_stx)+[xloc,yloc] color=255 plots,[llxy[0],ulxy[0]],[llxy[1],ulxy[1]],thick=thick,color=color,/device plots,[llxy[0],lrxy[0]],[llxy[1],lrxy[1]],thick=thick,color=color,/device plots,[ulxy[0],urxy[0]],[ulxy[1],urxy[1]],thick=thick,color=color,/device plots,[lrxy[0],urxy[0]],[lrxy[1],urxy[1]],thick=thick,color=color,/device ; get annotated image im=tvread(/greyscale,/quiet) ; optionally rotate to up=up for STX location if (rotate_up) then begin reformimage,im,imrot,rotate=angle_limb im=imrot endif ; optionally cut xrange and yrange if (xrange[1] ne -1 and yrange[1] ne -1) then begin if (min([[xrange],[yrange]]) lt 0 $ or max([[xrange],[yrange]]) gt 4095) then begin print,' ##### sdo_stx_fulldisk abort: wrong xrange or yrange' print,' xrange = '+trimd(xrange)+' yrange = '+trimd(yrange) retall endif im=im[xrange[0]:xrange[1],yrange[0]:yrange[1]] endif ; write image file writefits,resultdir+'/'+imagefile,im ; clean wdelete,22 end ; =============== main for testing per IDLWAVE H-c ====================== cd,'/home/rutten/data/SST/2016-09-05-demo' sdowav='171' stxname='sst' markcenter=1 marknorth=1 rotate_up=1 ;; xrange=[1500,2500] ;; yrange=[1500,2500] ; copy call at top sdo_stx_fulldisk,sdowav,stxname,$ rotate_up=rotate_up,markcenter=markcenter,marknorth=marknorth,$ xrange=xrange,yrange=yrange,thick=thick,$ alignfile=alignfile,$ fulldiskdir=fulldiskdir,level2dir=level2dir,resultdir=resultdir,$ datetime=datetime,solx=solx,soly=soly,nx_stx=nx_stx,ny_stx=ny_stx,$ px_stx=px_stx,angle_stx=angle_stx ; inspect showex,/allsdo,sdodir='fulldisk/annot' end