; file: sdo_diskfigs.pro = get full-disk SDO data and make figures ; init: Jun 24 2020 Rob Rutten Deil from wrt books cdej rob makefigssdo ; last: Dec 12 2020 Rob Rutten Deil ; note: can be run for every minute since SDO 2010 start until week ago ; todo: also scatcont with delays ;+ pro sdo_diskfigs,datetimes,sdowavs,gongha=gongha,$ wavpairs=wavpairs,outdir=outdir,newfiles=newfiles,ps2eps=ps2eps,$ ; imagefigs and options imagefigs=imagefigs,sdocolors=sdocolors,border=border,addaxes=addaxes,$ xbox=xbox,ybox=ybox,boxsizex=boxsizex,boxsizey=boxsizey,$ addlabel=addlabel,diskzooms=diskzooms,$ ; tilefigs and options tilefigs=tilefigs,diffheights=diffheights,$ arrowsize=arrowsize,plotzonal=plotzonal ; download SDO full-disk images and plot diverse figures from them ; ; INPUT PARAMETERS: ; datetimes = string array 1 or more 'YYYY.MM.DD_HH:MM:SS' strings ; sdowavs = string array 1 or more SDO wavnames as '304', 'magnetogram' ; ; OPTIONAL KEYWORD INPUTS: ; gongha = 1/0: also GONG Halpha (although pretty lousy) ; wavpairs: array of 2-elem string arrays for tile plots, either ; integer pairs = RR SDO wavnumbers, or string pairs = level2 names ; outdir: to write images+results in, default /tmp may overflow ; newfiles = 1/0: get new fits files even if already present (0) ; ps2eps = 1/0: convert ps output to eps for publication (0; takes long) ; imagefigs = 1/0: write SDO image (e)ps files in outdir (1) ; sdocolors = 1/0: use Schrijver SDO colors instead of monochrome (0) ; border: remove edge over this width, defaults 150 or 350 (below) ; addaxes = 1/0: add (X,Y) axes (0) ; xbox,ybox,boxsizex,boxsizey: add box for sdo_getdata_rr subfield ; addlabel = 1/O: add date-time-wav label in lower-left corner (0) ; diskzooms = 1/0: add magnified cutouts in corners (CdeJ style) (0) ; tilefigs = 1/0: write SDO tile ps files (e)ps in outdir (0) ; diffheights: formation height difference; default -1 = from muck ; arrowsize: magnifier to chart arrow length (default = adapt) ; plotzonal = 1/0: plot zonally-averaged shift sizes + radial component ; ; PROBLEMS: ; /tmp may overload for many datetimes, then set outdir to datadisk ; JSOC kills when busy: do NOT run simultaneous with sdo_getdata_rr.pro ; ; OUTPUTS: ; image fits files and (e)ps figure files in outdir or /tmp ; ; HISTORY: ; Jul 26 2020 RR: start ;- ; answer no-parameter query if (n_params(0) lt 2) then begin sp,sdo_diskfigs return endif ; set wall-clock timer (seconds) timestart=systime(1) ; keyword defaults if (n_elements(gongha) eq 0) then gongha=0 if (n_elements(outdir) eq 0) then outdir='/tmp' if (n_elements(newfiles) eq 0) then newfiles=0 if (n_elements(ps2eps) eq 0) then ps2eps=0 if (n_elements(imagefigs) eq 0) then imagefigs=1 ; !! whyse not if (n_elements(sdocolors) eq 0) then sdocolors=0 if (n_elements(border) eq 0) then border=150 ; !! sticky parameter if (n_elements(addaxes) eq 0) then addaxes=0 if (n_elements(xbox) eq 0) then xbox=0 if (n_elements(ybox) eq 0) then ybox=0 if (n_elements(boxsizex) eq 0) then boxsizex=0 if (n_elements(boxsizey) eq 0) then boxsizey=0 if (n_elements(addlabel) eq 0) then addlabel=0 if (n_elements(diskzooms) eq 0) then diskzooms=0 if (n_elements(tilefigs) eq 0) then tilefigs=0 if (n_elements(wavpairs) eq 0) then wavpairs='' if (n_elements(diffheights) eq 0) then diffheights=-1 ; special, maybe 0 if (n_elements(arrowsize) eq 0) then arrowsize=0 ; adaptive if (n_elements(plotzonal) eq 0) then plotzonal=0 ; check format of datetimes ndatetimes=n_elements(datetimes) for idatetime=0,ndatetimes-1 do begin if (strmid(datetimes[idatetime],4,1) ne '.' or $ strmid(datetimes[idatetime],7,1) ne '.' or $ strmid(datetimes[idatetime],10,1) ne '_' or $ strmid(datetimes[idatetime],13,1) ne ':') then begin print, ' ##### sdo_diskfigs abort: datetime '+$ datetimes[idatetime]+ ' not YYYY.MM.DD_HH:MM:SS' return endif endfor ; check axes versus zooms if (diskzooms eq 1 and addaxes eq 1) then begin print,' ===== sdo_diskfigs warning: no axes when zoomins' addaxes=0 endif ; diskzooms need specific border if (diskzooms eq 1 and addaxes eq 0) then border=350 ; set pixels per tile pxtile=50 ; im_patch box requested if (boxsizex gt 0 and boxsizey gt 0) then begin addbox=1 addaxes=1 endif else addbox=0 ; define permitted level2 fitsfilename strings (rrwavnr = index+1: ; 94,131,171,193,211,304,335,1600,1700,4500,mag,cont,dop ; 1 2 3 4 5 6 7 8 9 10 11 12 13 wavs=['94','131','171','193','211','304','335',$ ; rrwav 1-7 = AIA EUV '1600','1700','4500',$ ; rrwav 8-10 = AIA UV 'magnetogram','continuum','Dopplergram'] ; rrwav 11-13 = HMI nwavs=n_elements(wavs) ; check input sdowavs and fill rrwav array nsdowavs=n_elements(sdowavs) rrwav=intarr(nsdowavs) for isdowav=0,nsdowavs-1 do begin iwav=-1 for iw=0,nwavs-1 do if (wavs[iw] eq sdowavs[isdowav]) then iwav=iw if (iwav eq -1 ) then begin print,' ===== ERROR: sdowav invalid = ',sdowavs[isdowav] print,' ===== valid: ',wavs return endif rrwav[isdowav]=iwav endfor ; define Schrijver sdocolor wav numbers in RR order colorwavs=[94,131,171,193,211,304,335,1600,1700,1,1,1] ; start outer big loop over datetimes for idatetime=0,ndatetimes-1 do begin ; get GONG Halpha for this datetime if (gongha eq 1) then begin hafile='' if (newfiles eq 1) then hafile=gong_getimage(datetimes[idatetime],$ outdir=outdir) $ else begin if (file_test(outdir+'/gong_ha_'+datetimes[idatetime]+'*.fits') eq 1) $ then files=$ file_search(outdir+'/gong_ha_'+datetimes[idatetime]+'*.fits') $ else files=gong_getimage(datetimes[idatetime],outdir=outdir) hafile=files[0] endelse if (hafile[0] ne '' and hafile[0] ne -1) then $ print,' ----- sdo_diskfigs uses '+hafile[0] $ else begin print,' ##### no GONG Halpha for '+datetimes[idatetime]+'; skipped' goto, GONGEND endelse imha=readfits(hafile,hahead) gongtime=fxpar(hahead,'date_obs') ; eg: '2020-04-27T03:00:10' bytimha=bytscl(imha) ; AHA GONG keyword radius is always 900 px: all images scaled to that ; ?? later 860? ; and centered with [1024,1024] = disk center ; https://www.nso.edu/data/nisp-data/h-alpha/ datetimeut=anytim2utc(datetimes[idatetime],/ccsds) rsun_obs=get_rb0p(datetimeut,/radius,/quiet) ;; arcsecpx=rsun_obs/900. arcsecpx=rsun_obs/860. nx=2048 ny=2048 ; get (X,Y) axes if requested (no inserts) if (addaxes eq 1) then begin xcen=0 ycen=0 xoffset=nx/2 yoffset=ny/2 xaxisarr=(indgen(nx)*float(nx)/(nx-1)-xoffset)*arcsecpx+xcen yaxisarr=(indgen(ny)*float(ny)/(ny-1)-yoffset)*arcsecpx+ycen endif ; make ps plot psfilename=outdir+'/gong_ha_'+datetimes[idatetime]+'.ps' if (addaxes eq 0) then begin plotarea=[0.,0.,1.,1.] ; no white borders since no axes xtitle='' ytitle='' fontsize=9 endif else begin plotarea=[0.13,0.13,0.95,0.95] ; must be square xtitle='X [arcsec]' ytitle='Y [arcsec]' fontsize=5 ; maybe larger for multi-image figures endelse openpsplot,psfilename,thick=2,fontsize=fontsize,xsize=8.8,ysize=8.8 tv,bytimha,plotarea[0],plotarea[1],$ xsize=plotarea[2]-plotarea[0],ysize=plotarea[3]-plotarea[1],/normal ; add datetime label when requested if (addlabel eq 1) then $ xyouts,plotarea[0]+0.02,plotarea[1]+0.02,/norm,color=255,$ datetimes[idatetime]+' GONG Ha' ; add axes when requested if (addaxes eq 1) then contour,bytimha,xaxisarr,yaxisarr,$ /nodata,/noerase,/xstyle,/ystyle,position=plotarea,$ xticklen=-0.03,yticklen=-0.03,xtitle=xtitle,ytitle=ytitle ; add im_patch cutout box when specified via boxsizex and boxsizey if (addbox eq 1) then tvbox,[boxsizex,boxsizey],xbox,ybox,255,/data closepsplot,psfilename,opengv=0 if (ps2eps eq 1) then spawn,'pscropepsone '+psfilename print,' ----- wrote '+psfilename GONGEND: endif ; --------------- all done if only GONG requested if (imagefigs eq 0 and tilefigs eq 0) then return ; --------------- get SDO image files for this datetime ; define all-image array allimages=fltarr(4096,4096,nsdowavs) ; define all-time array alltimes=strarr(nsdowavs) ; SDO pixel size arcsecpx=0.6 ; AIA px size holds also for HMI after aia_prep below ; start inner big loop over sdowavs for isdowav=0,nsdowavs-1 do begin ; set psfilename filename=outdir+'/sdo_'+datetimes[idatetime]+'_'+sdowavs[isdowav] fitsfilename=filename+'.fits' psfilename=filename+'.ps' ; get image, new from JSOC or already gotten from outdir if (newfiles eq 1) then $ sdo_getimage,datetimes[idatetime],sdowavs[isdowav],im,indexim, $ outdir=outdir $ else begin if (file_test(fitsfilename) eq 1) then begin sdo_getimage,fitsfilename,sdowavs[isdowav],im,indexim,$ outdir=outdir print,' ----- sdo_diskfigs uses '+fitsfilename endif else begin sdo_getimage,datetimes[idatetime],sdowavs[isdowav],im,indexim,$ outdir=outdir endelse endelse ; keywords xcen=indexim.xcen ycen=indexim.ycen ; store observing time alltimes[isdowav]=indexim.t_obs ; pretty up if (rrwav[isdowav] lt 10) then im=histo_opt_rr(im,1E-4) if (rrwav[isdowav] eq 8 or rrwav[isdowav] eq 9) then im=sqrt(im>1) ;; ; check im ;; sv,sqrt(histo_opt_rr(im,1E-3)) ;; STOP ; enter this image in all-image array for current datetime allimages[*,*,isdowav]=im[*,*] ; ========== make image fig for this image, with or without zoom inserts if (imagefigs eq 1) then begin ; optional Schrijver SDO colors if (sdocolors eq 1 and isdowav lt 9) then $ aia_lct,r,g,b,wavelnth=colorwavs[rrwav[isdowav]],/load else $ aia_lct,r,g,b,wavelnth=colorwavs[9],/load ;; ; check ;; print,' ----- color: ',rrwav[isdowav],colorwavs[rrwav[isdowav]] ; crop borders im=im[border:4096-border,border:4096-border] sizeim=size(im) nx=sizeim[1] ny=sizeim[2] ; bytscl to make 255 "white" (mag maintain zero) bytim=bytscl(im) if (sdowavs[isdowav] eq 'magnetogram') then begin im=im>(-100)<(100) ; clips in "apparent" Gauss bytim=bytzero(im) endif ; ====== diskzooms set: add cuts with magnified inserts in corners if (diskzooms eq 1 and addaxes eq 0) then begin imins=bytim ; will be sucessively filled with cutouts ; set cut dimensions cutsize=300 mag=3 ; loop over cutouts and magnified inserts for icut=0,3 do begin ; define the cut locations = center down to limb (shifts with border) if (icut eq 0) then cutcenter=[nx/2,ny/2] if (icut eq 1) then cutcenter=[nx/2,ny/2-500] if (icut eq 2) then cutcenter=[nx/2,ny/2-1000] if (icut eq 3) then cutcenter=[nx/2,+cutsize/2+130] cutxmin=cutcenter[0]-cutsize/2 cutxmax=cutcenter[0]+cutsize/2 cutymin=cutcenter[1]-cutsize/2 cutymax=cutcenter[1]+cutsize/2 ; get cut image, bytscl independently to show fine structure best cutim=im[cutxmin:cutxmax,cutymin:cutymax] cutim=scale(cutim,mag,mag) ; ancient Herbie caller to rebin if (sdowavs[isdowav] ne 'magnetogram') then cutim=bytscl(cutim) $ else cutim=bytzero(cutim) ; add white borders to inserts thick=5 sizecut=size(cutim) borderim=bytarr(sizecut[1]+2*thick,sizecut[2]+2*thick)+255 ; all white sizeborder=size(borderim) nxins=sizeborder[1] nyins=sizeborder[2] borderim[thick:nxins-1-thick,thick:nyins-1-thick]=cutim[*,*] ; stick inserts in corners if (icut eq 0) then imins[0:nxins-1,ny-nyins:ny-1]=borderim[*,*] if (icut eq 1) then imins[nx-nxins:nx-1,ny-nyins:ny-1]=borderim[*,*] if (icut eq 2) then imins[nx-nxins:nx-1,0:nyins-1]=borderim[*,*] if (icut eq 3) then imins[0:nxins-1,0:nyins-1]=borderim[*,*] ; add white cutout frames thick=4 ; in px imins[cutxmin:cutxmax,cutymin-thick:cutymin+thick]=255 imins[cutxmin:cutxmax,cutymax-thick:cutymax+thick]=255 imins[cutxmin-thick:cutxmin+thick,cutymin:cutymax]=255 imins[cutxmax-thick:cutxmax+thick,cutymin:cutymax]=255 endfor ; end of loop over cutout inserts endif else imins=bytim ; no inserts ;; ; inspect ;; showex,imins ;; STOP ; get axes if requested (no inserts) if (addaxes eq 1) then begin xoffset=nx/2 yoffset=ny/2 xaxisarr=(indgen(nx)*float(nx)/(nx-1)-xoffset)*arcsecpx+xcen yaxisarr=(indgen(ny)*float(ny)/(ny-1)-yoffset)*arcsecpx+ycen endif ; make ps plot if (addaxes eq 0) then begin plotarea=[0.,0.,1.,1.] ; no white borders since no axes xtitle='' ytitle='' fontsize=9 endif else begin plotarea=[0.13,0.13,0.95,0.95] ; must be square xtitle='X [arcsec]' ytitle='Y [arcsec]' fontsize=5 ; maybe larger for multi-image figures endelse openpsplot,psfilename,thick=2,fontsize=fontsize,xsize=8.8,ysize=8.8 tv,imins,plotarea[0],plotarea[1],$ xsize=plotarea[2]-plotarea[0],ysize=plotarea[3]-plotarea[1],/normal ; add datetime label when requested if (addlabel eq 1) then $ xyouts,plotarea[0]+0.02,plotarea[1]+0.02,/norm,color=255,$ datetimes[idatetime]+' '+sdowavs[isdowav] ; add axes when requested if (addaxes eq 1) then contour,imins,xaxisarr,yaxisarr,$ /nodata,/noerase,/xstyle,/ystyle,position=plotarea,$ xticklen=-0.03,yticklen=-0.03,xtitle=xtitle,ytitle=ytitle ; add im_patch cutout box when specified via boxsizex and boxsizey if (addbox eq 1) then tvbox,[boxsizex,boxsizey],xbox,ybox,255,/data closepsplot,psfilename,opengv=0 if (ps2eps) then spawn,'pscropepsone '+psfilename print,' ----- wrote '+psfilename ; back to monochrome tvlct,byte(findgen(256)),byte(findgen(256)),byte(findgen(256)) endif ; end of imagefigs production for this image endfor ; end loop over sdowavs ; ===== make tileplots when requested and tilepair(s) given if (tilefigs ne 0) then begin ; check on 2 elements per pair szwavpairs=size(wavpairs) if (wavpairs[0] eq '' or szwavpairs[1] ne 2) then begin print,' ##### sdo_diskfigs: wrong wavpairs ='+trimd(wavpairs) STOP endif ; start loop over pairs npairs=n_elements(wavpairs)/2 for ipair=0,npairs-1 do begin wavpair=reform(wavpairs[*,ipair]) ; check wavpair szpair=size(wavpair) if (szpair[0] ne 1 or szpair[1] ne 2) then begin print,' ##### sdo_diskfigs: wavpair not 2-elem array:'+trimd(wavpair) STOP endif ; get rrpair values rrwavpair=[0,0] if (szpair[2] eq 2) then rrwavpair=wavpair ; integers = rrnumber pair if (szpair[2] eq 7) then begin ; strings = SDO level1 name pair for isdowav=0,nsdowavs-1 do begin if (sdowavs[isdowav] eq wavpair[0]) then $ rrwavpair[0]=rrwav[isdowav]+1 if (sdowavs[isdowav] eq wavpair[1]) then $ rrwavpair[1]=rrwav[isdowav]+1 endfor endif ; get images and times for isdowav=0,nsdowavs-1 do begin if sdowavs[isdowav] eq wavpair[0] then begin im1=allimages[*,*,isdowav] time1=alltimes[isdowav] endif if sdowavs[isdowav] eq wavpair[1] then begin im2=allimages[*,*,isdowav] time2=alltimes[isdowav] endif endfor ; ----- run imagepair mucking and tileshift finding ; muck this pair my way astile=pxtile*arcsecpx sdo_muckimagepair,im1,im2,rrwavpair,$ muckim1,muckim2,astile,tilecut,heightdiff ; use diffheight for heightdiff when specied in call if (diffheights[0] ne -1) then begin ndiffheights=n_elements(diffheights) if (ndiffheights ne npairs) then begin print,' ##### sdo_diskfigs: ndiffheights ne npairs' STOP endif heightdiff=diffheights[ipair] ; overwrite muckimagepair value endif ; set timediff for differential derotation (negligible anyhow) tim1tai=anytim2tai(time1) tim2tai=anytim2tai(time2) timediff=tim1tai-tim2tai ; tile-measure and write tileshifts shift=findimshift_tiled($ muckim1,muckim2,outdir=outdir,$ datetime=datetimes[idatetime],wavpair=wavpair,$ pxtile=pxtile,tilecut=tilecut,precision=0.001,$ heightdiff=heightdiff,timediff=timediff,$ writetileshifts=1,shiftchart=0,verbose=0,blink=0) ; make tilefigs maketilefigs,outdir=outdir,datetime=datetimes[idatetime],$ wavpair=wavpair,arcsecpx=arcsecpx,$ shiftchart=1,arrowsize=arrowsize,plotzonal=plotzonal endfor ; end loop over wavpairs endif ; end of tilefigs endfor ; end loop over datetimes ; print elapsed time if (ndatetimes gt 1) then begin timedone=systime(1) timelaps=ntostr((timedone-timestart)/60.,format='(F11.1)') print,' ===== sdo_diskfigs took '+timelaps+' minutes' endif if (outdir eq '/tmp') then $ print,' ===== wrote sdo_diskfigs results in /tmp: consider saving elsewhere' end ; end of program ; ============================================ run per Hyper C ;; ; sample all complete SDO years at start every month ;; outdir='/media/rutten/RRHOME/alldata/SDO/2020-UVtilecharts/2020-108samples' ;; years=['2011','2012','2013','2014','2015','2016','2017','2018','2019'] ;; months=['01','02','03','04','05','06','07','08','09','10','11','12'] ;; nyr=n_elements(years) ;; nmo=n_elements(months) ;; datetimes=strarr(nyr*nmo) ;; for iyr=0,nyr-1 do begin ;; for imo=0,nmo-1 do begin ;; datetimes[iyr*nmo+imo]=years[iyr]+'.'+months[imo]+'.01_00:00' ;; endfor ;; endfor ;; ; this took 249.3 min for 5 wavs, 1 pair, all figs => 74G, 2 GONGs missing datetimes='2020.01.01_00:00' ; for single-datetime tests ;; datetimes='2020.05.30_14:58:46' ; best-match SolO first 174 image ;; xbox=375+1 ; add 1 for rotation ;; ybox=148 ;; boxsizex=700 ;; boxsizey=700 ;; sdowavs=['magnetogram','1600','1700','304','171','193'] ;; sdowavs=['magnetogram','304','131'] ;; sdowavs=['1600','1700'] ; choices newfiles=1 gongha=0 imagefigs=0 sdocolors=0 addaxes=1 addlabel=1 diskzooms=0 tilefigs=1 wavpairs=['1600','1700'] ;; wavpairs=['304','magnetogram'] ;; wavpairs=['131','304'] diffheights=-1 ; use RR muck defaults diffheights=0 plotzonal=0 arrowsize=0 ; use defaults sdo_diskfigs,datetimes,sdowavs,gongha=gongha,$ outdir=outdir,newfiles=newfiles,ps2eps=ps2eps,$ imagefigs=imagefigs,$ ; imagefigs options sdocolors=sdocolors,border=border,addaxes=addaxes,$ xbox=xbox,ybox=ybox,boxsizex=boxsizex,boxsizey=boxsizey,$ addlabel=addlabel,diskzooms=diskzooms,$ tilefigs=tilefigs,$ ; tilefigs options wavpairs=wavpairs,diffheights=diffheights,$ arrowsize=arrowsize,plotzonal=plotzonal end ; immer 2 quadrupool zwarte plekken alsof quadrant niet deugt ; of fout px/arcsec schaling