;+ ; NAME: ; THM_ASI_CREATE_MOSAIC ; ; PURPOSE: ; create mosaic with all THEMIS ASI ; ; CATEGORY: ; None ; ; CALLING SEQUENCE: ; THM_ASI_CREATE_MOSAIC,time ; ; INPUTS: ; Time like '2006-01-01/05:00:00' ; ; OPTIONAL INPUTS: ; None ; ; KEYWORD PARAMETERS: ; cal_files calibration files if they do not need to be read ; pgm_file do not read CDF, but pgm-files ; verbose print some diagnostics ; insert insert stop before end of program ; ; gif_out create a gif-file ; gif_dir directory for gif-output ; ; exclude string of station names that should not be plotted ; show string of station names that should only be plotted ; minval minimum value for black ; maxval maximum value for white ; minimum_elevation minimum elevation to plot in degrees ; mask mask certain parts of image ; ; scale scale for map set ; central_lon geographic longitude of center of plot ; central_lat geographic latitude of center of plot ; rotation rotate map ; projection projection for map set, MAP_PROJ_INFO, PROJ_NAMES=names ; color_continent shade of continent fill ; color_background shade of background ; ; zbuffer do in z-buffer, not on the screen ; cursor finish with cursor info, loop if cursor>1 ; window set window number ; xsize xsize of window ; ysize ysize of window ; position=position position of plot on window (normal coordinates) ; noerase=noerase do not erase current window (no effect if {x,y}size set ; keep_z keep z-buffer open for further plotting ; ; no_grid=no_grid do not plot geomagnetic grid ; no_midnight=no_midnight do not plot midnight meridian ; no_label do not label mosaic with date and time ; add_plot stop because we want to add something ; force_map plot map even if there are no images ; ; xy_pos xy position ; location mark geographic location [lo,la] ; track1 mark geographic location [lo,la] ; track2 mark geographic location [lo,la] ; ; top top color to be used for polyfill ; no_color do not load color table, use existing ; xy_cursor create array of cursor selected values to pass to upper program ; ssize size of symbol for location ; sym_color color of location ; ; stoptime create multiple mosaics ; timestep time steps for multiple mosaics in seconds ; ; OUTPUTS: ; None ; ; OPTIONAL OUTPUTS: ; None ; ; COMMON BLOCKS: ; None ; ; SIDE EFFECTS: ; None ; ; RESTRICTIONS: ; None ; ; EXAMPLE: ; THM_ASI_CREATE_MOSAIC,'2006-01-01/05:00:00' ; THM_ASI_CREATE_MOSAIC,'2007-03-01/04:00:00',/thumb,exclude='ekat' ; ; MODIFICATION HISTORY: ; Written by: Harald Frey, 02/06/2007 ; based on example from Donovan/Jackel ; ; 2007-03-15, hfrey, thumbnails, keyword exclude ; 2007-03-27, hfrey, special treatment for not mapped thumbsnails ; 2007-12-21, jmm, added explicit set_plot,'z' for zbuffer ; 2008-07-21, jmm, added gif_dir, for output directory option ; 2009-06-17, hfrey, a few additions to make my life easier ; 2009-11-10, cgabrielse, added xy_cursor keyword for ; sending cursor values up level ; 2012-07-02, jmm, Added color_annotation keyword ; ; NOTES: ; ; VERSION: ; $LastChangedBy: hfrey $ ; $LastChangedDate: 2014-01-07 12:17:57 -0800 (Tue, 07 Jan 2014) $ ; $LastChangedRevision: 13811 $ ; $URL: svn+ssh://thmsvn@ambrosia.ssl.berkeley.edu/repos/spdsoft/trunk/projects/themis/ground/thm_asi_create_mosaic.pro $ ; ;- PRO THM_ASI_CREATE_MOSAIC,time,$ cal_files=cal_files,$ ; calibration files already read gif_out=gif_out,$ ; output in gif file verbose=verbose,$ ; print debug messages pgm_file=pgm_file,$ ; read raw pgm files thumb=thumb,$ ; force thumbnails over full images exclude=exclude,$ ; exclude certain stations top=top,$ ; set top value for image special=special,$ ; treat 32x32 thumbnails show=show,$ ; limit stations shown scale=scale,$ ; scale for map set central_lat=central_lat,$ ; geographic latitude of center of plot central_lon=central_lon,$ ; geographic longitude of center of plot color_continent=color_continent,$ ; shade of continent fill color_background=color_background,$; shade of background position=position,$ ; position of plot on window (normal coordinates) xsize=xsize,$ ; xsize of window ysize=ysize,$ ; ysize of window noerase=noerase,$ ; do not erase current window (no effect if {x,y}size set) zbuffer=zbuffer,$ ; do it in z-buffer cursor=cursor,$ ; finish with cursor info projection=projection,$ ; map projection maxval=maxval,$ ; brightness scaling of images minval=minval,$ ; brightness scaling of images window=window,$ ; set window number rotation=rotation,$ ; rotate map away from North up full_minute=full_minute,$ ; process a full minute of pgm files minimum_elevation=minimum_elevation,$ ; set minimum elevation merge=merge,$ ; merge full resolution and thumbnail images gif_dir=gif_dir,$ ; An output directory for the gif output, default is the local working dir. force_map=force_map,$ ; force display of empty map no_grid=no_grid,$ ; do not plot grid no_midnight=no_midnight,$ ; do not plot midnight meridian add_plot=add_plot,$ ; stop so we can add something mask=mask,$ ; mask part of image xy_pos=xy_pos,$ ; mark specific x,y-location location=location,$ ; mark geographic location=[long,lat] or [[lo,la],[lo,la],[lo,la]] no_color=no_color, $ ; do not load color table xy_cursor=xy_cursor,$ ; create an array that records the cursor output and passes it to upper program track1=track1,$ ; mark geographic location=[long,lat] or [[lo,la],[lo,la],[lo,la]] track2=track2,$ ; mark geographic location=[long,lat] or [[lo,la],[lo,la],[lo,la]] ssize=ssize,$ ; size of symbol for location sym_color=sym_color,$ ; color of location keep_z=keep_z,$ ; keep z-buffer open stoptime=stoptime,$ ; create multiple mosaics timestep=timestep,$ ; time step in seconds for multiple mosaics color_annotation = color_annotation,$ ; color for annotation, via xyouts insert=insert,$ ; insert stop before end of program no_label=no_label ; do not label mosaic with date and time ; input check if keyword_set(verbose) then verbose=systime(1) if (strlen(time) ne 19) then begin dprint, 'Wrong time input' dprint, 'Correct value like 2006-01-01/05:00:00' return endif if keyword_set(xy_pos) then begin dd=size(xy_pos) if (dd[2] ne 2) then begin dprint, 'XY_pos input wrong!' dprint, 'Needs to be given like [[x1,x2,x3,x4,x5,...],[y1,y2,y3,y4,y5,...]]' return endif endif ; check brightness scaling if keyword_set(maxval) then begin if not keyword_set(minval) then begin dprint, 'minval has to be set with maxval' return endif if not keyword_set(show) then begin dprint, 'Show has to be set with maxval' return endif if n_elements(show) ne n_elements(maxval) then begin dprint, 'N_elements of show and maxval have to match' return endif endif ; strip time res=time_struct(time) year=res.year month=res.month day=res.date hour=res.hour minute=res.min second=res.sec ; setup thm_init timespan,time,1,/hour thm_asi_stations,site,loc if keyword_set(zbuffer) then map_scale=2.6e7 else map_scale=4.e7 if keyword_set(scale) then map_scale=scale if not keyword_set(central_lon) then central_lon=255. if not keyword_set(central_lat) then central_lat=63. if not keyword_set(xsize) then xsize=700 if not keyword_set(ysize) then ysize=410 if not keyword_set(top) then top=254 ; characters if keyword_set(zbuffer) then chars=1.15 else chars=1.5 ; some setup if keyword_set(minimum_elevation) then minimum_elevation_to_plot=minimum_elevation else minimum_elevation_to_plot=8. ;degrees if (keyword_set(thumb) and not keyword_set(special)) then n1=1024l else n1=256l*256l ;color for annotation If(keyword_set(color_annotation)) Then Begin xyouts_color = color_annotation[0] Endif Else xyouts_color = 0 ; clean up before start names=tnames('thg_as*') if (names[0] ne '') then store_data,delete=names ; load cal_files once for the whole loop if keyword_set(stoptime) then begin if keyword_set(show) then thm_load_asi_cal,show,cal_files else $ thm_load_asi_cal,'atha chbg ekat fsmi fsim fykn gako gbay gill inuv kapu kian kuuj mcgr pgeo pina rank snkq tpas whit yknf nrsq snap talo',cal_files endif ; point to run another mosaic repeat_loop: ; search for midnight file if not keyword_set(no_midnight) then begin f=file_search('midnight.sav',count=midnight_count) if (midnight_count eq 1) then begin midlons=fltarr(40)+!values.f_nan ut_hour=float(hour)+float(minute)/60.+float(second)/3600. restore,f for i=0,39 do begin lon=interpol(findgen(141)+start_longitude,reform(midnight[i,*]),ut_hour) midlons[i]=lon[0] endfor bad=where(midlons gt 360.,count) if (count gt 0) then midlons[bad]=!values.f_nan endif endif else midnight_count=0 ; read available data thm_mosaic_array,year,month,day,hour,minute,second,strlowcase(site),$ image,corners,elevation,pixel_illuminated,n_sites,verbose=verbose,$ cal_files=cal_files,pgm_file=pgm_file,thumb=thumb,special=special,$ show=show,exclude=exclude,full_minute=full_minute,merge=merge,$ mask=mask ; exclude unwanted sites if keyword_set(exclude) then begin for i=0,n_elements(exclude)-1 do begin not_site=where(strlowcase(site) eq strlowcase(exclude[i]),count) if (count eq 0) then dprint, 'Not a valid site: ',exclude[i] else begin corners[*,*,*,not_site]=!values.f_nan image[*,not_site]=!values.f_nan endelse endfor endif ; fill variables bytval=fltarr(n_sites)+1. bitval=fltarr(n_sites) if keyword_set(maxval) then begin for i=0,n_elements(maxval)-1 do begin index=where(strlowcase(site) eq strlowcase(show[i])) bytval[index]=maxval[i] bitval[index]=minval[i] endfor for i=0,n_sites-1 do image[*,i]=bytscl(image[*,i],min=bitval[i],max=bytval[i]) endif else begin if keyword_set(special) then for i=0,n_sites-1 do bytval[i]=(median(image[0:1023,i]) > 1) else begin if keyword_set(full_minute) then begin for i=0,n_sites-1 do bytval[i]=(median(image[*,i,*]) > 1) for i=0,n_sites-1 do image[*,i,*]=((image[*,i,*]/bytval[i])*64.) < 254 endif else begin ; prevent divide by zero for i=0,n_sites-1 do bytval[i]=(median(image[*,i]) > 1) ; prevent divide by zero mm_dumm=where(merge ne -1,count_mm_dumm) if (count_mm_dumm ne 0) then begin for dki=0,count_mm_dumm-1 do bytval[mm_dumm[dki]]=(median(image[0:1023,mm_dumm[dki]]) > 1) endif for i=0,n_sites-1 do image[*,i]=((image[*,i]/bytval[i])*64.) < 254 endelse endelse endelse ; no images found if (max(bytval) eq 1.) then begin dprint, 'No images for ',time if not keyword_set(force_map) then begin if keyword_set(gif_out) then gif_out='' if double(!version.release) lt 8.0d then heap_gc return endif endif ; exclude unwanted sites if keyword_set(exclude) then begin for i=0,n_elements(exclude)-1 do begin not_site=where(strlowcase(site) eq strlowcase(exclude[i]),count) if (count eq 0) then dprint, 'Not a valid site: ',exclude[i] else begin corners[*,*,*,not_site]=!values.f_nan bytval[not_site]=!values.f_nan endelse endfor endif ;========================================================================= ; generate images for full minute if keyword_set(full_minute) then begin time_start=time ; run through images for ikk=0,19 do begin ;zbuffer needs to be set before the loadct call in thm_map_set, ;otherwise this bombs the second time through because of reset to 'x' ;later in this program, jmm 21-dec-2007 if(keyword_set(zbuffer)) then set_plot, 'z' ; set up the map thm_map_set,scale=map_scale,$ central_lat=central_lat,$ ; geographic latitude of center of plot central_lon=central_lon,$ ; geographic longitude of center of plot color_continent=color_continent,$ ; shade of continent fill color_background=color_background,$ ; shade of background position=position,$ ; position of plot on window (normal coordinates) xsize=xsize,$ ; xsize of window ysize=ysize,$ ; ysize of window noerase=noerase,$ ; do not erase current window (no effect if {x,y}size set zbuffer=zbuffer,$ projection=projection,$ window=window,$ rotation=rotation,$ no_color=no_color ; normal filling for pixel=0l,n1-1l do begin for i_site=0,n_sites-1 do begin if ((pixel_illuminated[pixel,i_site] eq 1) and $ (elevation[pixel,i_site] gt minimum_elevation_to_plot)) then $ polyfill,corners[pixel,[0,1,2,3,0],0,i_site],$ corners[pixel,[0,1,2,3,0],1,i_site],color=image[pixel,i_site, ikk] < top if keyword_set(special) then if ((spec_illu[pixel*64 < (n1-1),i_site] eq 1) and $ (elevation[pixel*64 < (n1-1),i_site] gt minimum_elevation_to_plot)) then begin for ijk=0,63 do polyfill,corners[pixel*64+ijk < (n1-1),[0,1,2,3,0],0,i_site],$ corners[pixel*64+ijk < (n1-1),[0,1,2,3,0],1,i_site],color=image[pixel*64+ijk < (n1-1),i_site] < top ; stop endif endfor endfor ; finish map return_lons=1 return_lats=1 if not keyword_set(no_grid) then thm_map_add,invariant_lats=findgen(4)*10.+50.,invariant_color=210,$ invariant_linestyle=1,/invariant_lons,return_lons=return_lons,$ return_lats=return_lats time=time_string(time_double(time_start)+ikk*3) if not keyword_set(no_label) then begin xyouts,0.005,0.018,time,color=xyouts_color,/normal,charsize=chars xyouts,0.005,0.060,'THEMIS-GBO ASI',color=xyouts_color,/normal,charsize=chars endif ; plot midnight file if (not keyword_set(no_midnight) and midnight_count eq 1) then $ plots,smooth(midlons-360.,5),findgen(40)+40.,color=255,/data ; gif output if keyword_set(gif_out) then begin If(keyword_set(gif_dir)) Then gdir = gif_dir Else gdir = './';jmm, 21-jul-2008 tvlct,r,g,b,/get img=tvrd() ; strip time because it changes through the loop res=time_struct(time) year=res.year month=res.month day=res.date hour=res.hour minute=res.min second=res.sec out_name='MOSA.'+year+'.'+month+'.'+day+'.'+hour+'.'+minute+'.'+second+'.gif' write_gif,gdir+out_name,img,r,g,b dprint, 'Output in ',out_name gif_out=out_name endif if keyword_set(zbuffer) then zbuffer=tvrd() if not keyword_set(keep_z) and keyword_set(zbuffer) then begin device,/close set_plot,'x' endif endfor endif else begin ; end of full_minute loop ;========================================================================= ;zbuffer needs to be set before the loadct call in thm_map_set, ;otherwise this bombs the second time through because of reset to 'x' ;later in this program, jmm 21-dec-2007 if(keyword_set(zbuffer)) then set_plot, 'z' ; set up the map thm_map_set,scale=map_scale,$ central_lat=central_lat,$ ; geographic latitude of center of plot central_lon=central_lon,$ ; geographic longitude of center of plot color_continent=color_continent,$ ; shade of continent fill color_background=color_background,$ ; shade of background position=position,$ ; position of plot on window (normal coordinates) xsize=xsize,$ ; xsize of window ysize=ysize,$ ; ysize of window noerase=noerase,$ ; do not erase current window (no effect if {x,y}size set zbuffer=zbuffer,$ projection=projection,$ window=window,$ rotation=rotation,$ no_color=no_color ; fill it if keyword_set(special) then begin spec_illu=pixel_illuminated-2 for i=0,n_elements(special)-1 do begin spec_stat=where(strlowcase(site) eq strlowcase(special[i])) spec_illu[*,spec_stat]=pixel_illuminated[*,spec_stat] endfor endif ; normal filling for pixel=0l,n1-1l do begin for i_site=0,n_sites-1 do begin if ((pixel_illuminated[pixel,i_site] eq 1) and $ (elevation[pixel,i_site] gt minimum_elevation_to_plot)) then $ polyfill,corners[pixel,[0,1,2,3,0],0,i_site],$ corners[pixel,[0,1,2,3,0],1,i_site],color=image[pixel,i_site] < top if keyword_set(special) then if ((spec_illu[pixel*64 < (n1-1),i_site] eq 1) and $ (elevation[pixel*64 < (n1-1),i_site] gt minimum_elevation_to_plot)) then begin for ijk=0,63 do polyfill,corners[pixel*64+ijk < (n1-1),[0,1,2,3,0],0,i_site],$ corners[pixel*64+ijk < (n1-1),[0,1,2,3,0],1,i_site],color=image[pixel*64+ijk < (n1-1),i_site] < top endif endfor if (keyword_set(special) and pixel ge 1024) then break endfor ; finish map return_lons=1 return_lats=1 thm_map_add,invariant_lats=findgen(4)*10.+50.,invariant_color=210,$ invariant_linestyle=1,/invariant_lons,return_lons=return_lons,$ return_lats=return_lats,no_grid=no_grid if not keyword_set(no_label) then begin xyouts,0.005,0.018,time,color=xyouts_color,/normal,charsize=chars xyouts,0.005,0.060,'THEMIS-GBO ASI',color=xyouts_color,/normal,charsize=chars endif if keyword_set(verbose) then dprint, 'After map: ',systime(1)-verbose,$ ' Seconds' ; search for midnight file if (not keyword_set(no_midnight) and midnight_count eq 1) then $ plots,smooth(midlons-360.,5),findgen(40)+40.,color=255,/data ; stop so we can add something if keyword_set(add_plot) then stop ; mark ground tracks of satellites if keyword_set(track1) then begin plots,track1[0,*],track1[1,*],psym=3 endif if keyword_set(track2) then begin plots,track2[0,*],track2[1,*],psym=4 endif if keyword_set(location) then begin if keyword_set(ssize) then ssize=ssize else ssize=1 if keyword_set(sym_color) then scolor=sym_color else scolor=255 plots,location[0,*],location[1,*],psym=2,symsize=ssize,color=scolor endif if keyword_set(cursor) then begin ss=size(cursor) xy_cursor=fltarr(cursor,4) if (ss[1] ne 2) then cursor=1 for loop=1,cursor do begin dprint, 'Point cursor on map!' cursor,x,y,/data wait,0.25 res=convert_coord(x,y,/data,/to_device) dprint, 'Location: ',res,x,y xy_cursor[loop-1L,*]=[res[0],res[1],x,y] endfor endif ; input like [[x1,x2,x3,x4,x5,...],[y1,y2,y3,y4,y5,...]] if keyword_set(xy_pos) then begin dd=size(xy_pos) if (dd[0] eq 1) then begin res=convert_coord(xy_pos[0],xy_pos[1],/to_data,/device) dprint, 'Location: ',xy_pos,res[0:1],format='(a12,2i5,2f10.3)' xy_pos_out=[xy_pos[0],xy_pos[1],res[0],res[1]] endif else begin xy_pos_out=fltarr(dd[1],4) res=convert_coord(xy_pos[*,0],xy_pos[*,1],/to_data,/device) for i1=0L,dd[1]-1L do begin dprint, 'Location: ',xy_pos[i1,*],res[0:1,i1],format='(a12,2i5,2f10.3)' xy_pos_out[i1,*]=[xy_pos[i1,0],xy_pos[i1,1],res[0,i1],res[1,i1]] endfor endelse xy_pos=xy_pos_out endif ; gif output if keyword_set(gif_out) then begin If(keyword_set(gif_dir)) Then gdir = gif_dir Else gdir = './' tvlct,r,g,b,/get img=tvrd() ; now add the secret code of input parameters img[40:43,0]=[13,251,117,239] ; time of mosaic img[0:6,0]=[year/100,year-(year/100)*100,month,day,hour,minute,second] ; thumb flag if keyword_set(thumb) then img[7,0]=1 else img[7,0]=0 ; central_lon and lat of mosaic if (central_lon lt 0.) then central_lon=central_lon+360. img[8:12,0]=[fix(central_lon)/100,fix(central_lon)-(fix(central_lon)/100)*100,$ fix((central_lon-fix(central_lon))*100),$ fix(central_lat),fix((central_lat-fix(central_lat))*100)] ; map_scale res=strsplit(string(map_scale*1.e10),'e',/extract) img[13:15,0]=[fix(res[0]),fix((float(res[0])-fix(res[0]))*100),fix(res[1])-10] ; xsize and ysize img[16:19,0]=[xsize/100,xsize-(xsize/100)*100,ysize/100,ysize-(ysize/100)*100] ; rotation if keyword_set(rotation) then begin if (rotation lt 0.) then rotation=rotation+360. img[20:22,0]=[fix(rotation/100),fix(rotation-(fix(rotation)/100)*100),$ fix((rotation-fix(rotation))*100)] endif else img[20:22]=[0,0,0] ; minimum elevation img[23:24,0]=[fix(minimum_elevation_to_plot),$ fix((minimum_elevation_to_plot-fix(minimum_elevation_to_plot))*100)] ; zbuffer if keyword_set(zbuffer) then img[25,0]=1 else img[25,0]=0 ; code stations img[49,0]=n_sites for i1=0,n_sites-1 do begin case 1 of finite(bytval[i1]) eq 0: img[50+i1,0]=0 bytval[i1] eq 1.: img[50+i1,0]=1 bytval[i1] gt 1.: img[50+i1,0]=2 endcase endfor ; construct the name out_name='MOSA.'+string(year,'(i4.4)')+'.'+string(month,'(i2.2)')+'.'+$ string(day,'(i2.2)')+'.'+string(hour,'(i2.2)')+'.'+string(minute,'(i2.2)')+$ '.'+string(second,'(i2.2)')+'.gif' write_gif,gdir+out_name,img,r,g,b dprint, 'Output in ',out_name gif_out=out_name tv,img endif if keyword_set(zbuffer) then zbuffer=tvrd() if not keyword_set(keep_z) and keyword_set(zbuffer) then begin device,/close set_plot,'x' endif endelse ; single time ;========================================================================= ; loop of mosaics if keyword_set(stoptime) then begin if keyword_set(timestep) then new_time=time_double(time)+timestep else new_time=time_double(time)+3.d0 ; strip time res=time_struct(new_time) year=res.year month=res.month day=res.date hour=res.hour minute=res.min second=res.sec if (new_time le time_double(stoptime)) then begin time=time_string(new_time) goto,repeat_loop endif endif if keyword_set(verbose) then dprint, 'Calculation took ',systime(1)-verbose,$ ' Seconds' if double(!version.release) lt 8.0d then heap_gc if keyword_set(insert) then stop end