;+ ; PROJECT: Automated Region Selection Extraction ; alternative (J.Ireland) Forward Extraction Collection of Knowledge ; ; FUNCTION: arse2struct ; ; PURPOSE: Read in MDI full disc map and return ARSE parameters ; ; USEAGE: result = arse2struct( map ) ; ; INPUT: ; data MDI full disc map ; ; KEYWORDS: NONE ; ; OUTPUT: ; result - A structure containing ARSE paramters list ; ; This is a structure with one entry for each selected region ; The entries are ; 'Time_obs' - time of data ; 'ARSE_data' - projection corrected data ; 'ARSE_contour' - contoured region. ; Each filled contour has a different value [0,-10,...,-10*n] ; So it is possible to select different contours ; 'ARSE_nonpixels' - 1 where pixels belong to other active regions, else 0 ; also if there are a full row/column of 1, there is probably dropped data ; 'ARSE_num' - region number from ARSE, in decreasing order of maximum intensity pixel ; 'ARSE_posn_EACP', - position on projected full disc data ; 'ARSE_Size_Units', - Units of size, Mm^2 ; 'ARSE_Pos_Size', - size of all positive field ; 'ARSE_Neg_Size', - size of all negative field ; 'ARSE_Total_Size', - total size ; 'ARSE_Flux_Units', - units of flux, 10^15Mx' ; 'ARSE_Pos_Flux', - flux from positive field ; 'ARSE_Neg_Flux', - flux from negative field ; 'ARSE_Abs_Flux', - total absolute flux ; 'ARSE_Total_Flux', - total flux, with +ve and -ve cancelling, flux inbalance ; 'ARSE_Mixing', - [number of mixed contours, total number of contours] ; 'NOAA_AR', - NOAA designated number, if any ; 'NOAA_Size', - NOAA measured size, if any ; 'NOAA_Class', - NOAA Mount Wilson class, if any ; 'NOAA_C_flares' - list of C-flares from this region on this day, if any ; 'NOAA_C', - total number of C-class flares from this region on this day ; 'NOAA_M_flares',- list of M-flares from this region on this day, if any ; 'NOAA_M', - total number of M-class flares from this region on this day ; 'NOAA_X_flares',- list of X-flares from this region on this day, if any ; 'NOAA_X', - total number of X-class flares from this region on this day ; 'NOAA_Posn_HG' - NOAA heliographic position ; 'NOAA_Posn_HC' - NOAA heliocentric position ; 'NOAA_Posn_EACP' - NOAA projected corrected position ; 'Message', - Pertinent messages relating to ARSE extraction and/or NOAA region and/or flares. ; ; ;NB: IF AR selected has very low values it is may be an artificial spike in the image (e.g. cosmic ray hit) ; ; ; ; EXAMPLE: arse_struct = arse2struct( 'mdi_map' ) ; ; AUTHOR: 17-Feb-2005 R.T.James McAteer ; VERSION 1.1, changed so that if the data fails contouring, or region-growing ; it removes the window around this pixel and seeks a new seed pixel. previously it failed ; contouring or region-growing I declared end of algorithm. However I want to ; include data closer to the limb (h_angle=80), hence there will be more streaks ; in the image. ; ALSO, if the seed pixel is on the outer 10% then window around pixel is removed, and it ; seeks a new seed pixel. ; ALSO introduce a minimum field strength of 500G. If there are no pixels left below this value ; then the algorithm finished ;KNOWN BUGS (1)sometimes this will identify 2 close regions as 2 ARSE region. Only one NOAA number will be ; assigned to this ARSE region, which may cause errors (if, fo instance, the unselected region flared) ; (2) with the extended h_angle, it is possible that a region very close to the edge (but not in the outer ; 10% will be grew. hence be very careful about counting these as real regions. Even for real regions this ; edge data is going be strange (large fields, due to large cosine, and lots of streaks). Maybe I should ; just remove this completely? could just set the outer 10% to zero, or flag the >h_angle pixels in ; mdi_cos_correct and then set these flagged pixels to zero. ;- FUNCTION arse2struct, map, display=display ;correct for cos & mercator projection map=mdi_cos_correct(map,h_angle=80) ;get flare and NOAA stats flare_stats,map,xypos,noaa,flares,err=err_flare_stats ;carry out the cyclindrical equal area projection carmap=map_carrington3(map.data,[map.xc,map.yc],$ map.time,map.dx,xypos) sz_x=(size(carmap))(1) sz_y=(size(carmap))(2) ;make a copy of the carmap copymap=carmap ;make a blank array for contours contmap=fltarr(sz_x,sz_y) ;make a blank array for regions regmap=fltarr(sz_x,sz_y) ;select a window size ;this parameter is very important ;if it is too small, there will be lots of vertical/horizontal lines in the arse.data (i.e. the AR ;extends outside the window). if is it too large it is difficult to seperate regions close to each other ;so note that this number created an artifical maximum area of (win_size)^2 win_size=150 ;make a structure to contain the ARSE results arse_struct = create_struct( 'Time_obs','',$ 'ARSE_data', FLTARR(2*win_size,2*win_size),$ 'ARSE_contour', FLTARR(2*win_size,2*win_size),$ 'ARSE_nonpixels', FLTARR(2*win_size,2*win_size),$ 'ARSE_num', -1,$ 'ARSE_posn_EACP', [0.,0.],$ 'ARSE_Size_Units','Mm^2',$ 'ARSE_Pos_Size', 0.,$ 'ARSE_Neg_Size', 0.,$ 'ARSE_Total_Size', 0.,$ 'ARSE_Flux_Units','10^15Mx',$ 'ARSE_Pos_Flux', 0.,$ 'ARSE_Neg_Flux', 0.,$ 'ARSE_Abs_Flux', 0.,$ 'ARSE_Total_Flux', 0.,$ 'ARSE_Mixing',[0.,0.],$ 'NOAA_AR','-1',$ 'NOAA_Size',0.,$ 'NOAA_Class','-1.',$ 'NOAA_Mac','-1.',$ 'NOAA_A_flares',FLTARR(15),$ ;allows for 15 flares 'NOAA_A',0,$ 'NOAA_B_flares',FLTARR(15),$ ;allows for 15 flares 'NOAA_B',0,$ 'NOAA_C_flares',FLTARR(15),$ ;allows for 15 flares 'NOAA_C',0,$ 'NOAA_M_flares',FLTARR(15),$ ;allows for 15 flares 'NOAA_M',0,$ 'NOAA_X_flares',FLTARR(15),$ ;allows for 15 flares 'NOAA_X',0,$ 'NOAA_Posn_HG',INTARR(2),$ 'NOAA_Posn_HC',FLTARR(2),$ 'NOAA_Posn_EACP',FLTARR(2),$ 'Message','') output = replicate(arse_struct, 50) ;number the first active region ar_number=0 already_identified = ['1'] ;find the next active region find_next_region: ;set the time of observation output[ar_number].Time_obs = map.time ;changed 090305, jma ;if the region cannot be contoured, or grew, or is too close to the edge, then go here cannot_contour_image: no_pixels_in_contour: skip_printing_stats: ;this chosen cutoff value is arbitrary (as is the chosen total number of selected regions, 50). ;however a cut off is needed to keep the program running quickly. IF max(abs(carmap)) LT 500 THEN GOTO, no_regions_remaining ;select AR seed pixel imagemax=where(abs(carmap) eq max(abs(carmap))) imagemax=imagemax(0) printxy, carmap, imagemax, x_p, y_p ;test code added 090305, jma IF KEYWORD_SET(display) THEN BEGIN print,ar_number,max(abs(carmap)) loadct,0 tvim,carmap,range=[-500,500] loadct,3 mark,x_p,y_p,psym=7,color=128 ENDIF ;need to extract a reasonable guess submap submap=carmap(0>(x_p-win_size):(sz_x-1)<(x_p+win_size-1),0>(y_p-win_size):(sz_y-1)<(y_p+win_size-1)) subimagemax=where(abs(submap) eq max(abs(submap))) subimagemax=subimagemax(0) printxy, submap, subimagemax,xs_p,ys_p ;IF AR selected is too close to egde of image. it is probably a projection correction effect. IF [x_p LT [fix(0.1*sz_x)] OR x_p GT fix(0.90*sz_x) OR y_p LT fix(0.1*sz_y) OR y_p GT fix(0.9*sz_y)] THEN BEGIN carmap(0>(x_p-win_size):(sz_x-1)<(x_p+win_size-1),0>(y_p-win_size):(sz_y-1)<(y_p+win_size-1))=0 ;output[ar_number].Message=[output[ar_number].Message + 'Too close to egde '] GOTO, skip_printing_stats ENDIF ;contour the region using 'contig_pixels' cor_cont=contig_pixels(submap,thresh=50,/choose_params2,centxy=[xs_p,ys_p],/nodisp,/quiet) IF max(cor_cont) eq 0 THEN BEGIN carmap(0>(x_p-win_size):(sz_x-1)<(x_p+win_size-1),0>(y_p-win_size):(sz_y-1)<(y_p+win_size-1))=0 ;output[ar_number].Message=[output[ar_number].Message + 'Cannot contour image '] GOTO, cannot_contour_image ENDIF ;so now cor_cont is the contours of the AR in the submap ;grow the submap region using 'calc_reg_growing3' and cor_cont as the boundary calc_reg_growing3,submap,cor_cont,filled,filled_contour,/quiet IF min(filled) GT -10 THEN BEGIN carmap(0>(x_p-win_size):(sz_x-1)<(x_p+win_size-1),0>(y_p-win_size):(sz_y-1)<(y_p+win_size-1))=0 ;output[ar_number].Message=[output[ar_number].Message + 'Cannot fill contour '] GOTO, no_pixels_in_contour ENDIF IF err_flare_stats NE '' THEN BEGIN output[ar_number].Message = [output[ar_number].Message] + err_flare_stats GOTO,cannot_get_NAR_info ENDIF ;check to see if any AR exists ;IF err_flare_stats eq 'NOAA data not found for specified times' THEN BEGIN ; output[ar_number].Message = [output[ar_number].Message] + err_flare_stats ; GOTO,cannot_get_NAR_info ;END ;check to see if any AR exist ;IF err_flare_stats eq 'All points rotated off disk' THEN BEGIN ; output[ar_number].Message = [output[ar_number].Message] + err_flare_stats ; GOTO,cannot_get_NAR_info ;END dist = REFORM(sqrt( (x_p-xypos(0,*))^2 + (y_p-xypos(1,*))^2)) noaa_identified = where(dist eq MIN(dist)) ;FOR i=0, N_ELEMENTS(dist)-1 DO BEGIN ; IF dist(i) LT 50 THEN BEGIN ; noaa_identified = i ; GOTO, got_region ; ENDIF ;ENDFOR ; ;got_region: ; IF min(dist) GT 100. THEN BEGIN noaa_reg = -1. output[ar_number].Message = [output[ar_number].Message +'Cannot identify NOAA region. '] GOTO, cannot_get_NAR_info ENDIF ELSE BEGIN noaa_reg_int = MIN(noaa (noaa_identified)) ; to account for duplication of NOAA regions noaa_identified = MIN(noaa_identified) ;this was a get_nar problem, which I may have solved in get_nar2 ENDELSE IF noaa_reg_int LE 9 THEN noaa_reg = '000'+ ARR2STR(noaa_reg_int,/trim) IF noaa_reg_int GE 10 AND noaa_reg_int LE 99 THEN noaa_reg = '00'+ ARR2STR(noaa_reg_int,/trim) IF noaa_reg_int GE 100 AND noaa_reg_int LE 999 THEN noaa_reg = '0'+ ARR2STR(noaa_reg_int,/trim) IF noaa_reg_int GE 1000 AND noaa_reg_int LE 9999 THEN noaa_reg = ARR2STR(noaa_reg_int,/trim) ;set the NOAA region number, if any IF (WHERE(already_identified eq noaa_reg))(0) ne -1 THEN GOTO, cannot_get_NAR_info output[ar_number].NOAA_AR = noaa_reg already_identified=[already_identified,noaa_reg] ;check to see if NOAA AR is identified with any flares ;********EXTENSION*********** ;NOTE I can add to flare_stats.pro and just identify flare positions ;so even if an AR is not numbered I can still identify flares from it ;**************************** IF err_flare_stats NE '' THEN BEGIN output[ar_number].Message=[output[ar_number].Message + err_flare_stats] GOTO, no_flares ENDIF flare_identified = where (flares.region EQ noaa_reg,count) IF count NE 0 THEN BEGIN classes = flares(flare_identified).fclass fluxes = flares(flare_identified).fbase a_flares = where (classes eq 'A',count) output[ar_number].NOAA_A = count IF count NE 0 THEN output[ar_number].NOAA_A_flares = fluxes(a_flares) b_flares = where (classes eq 'B',count) output[ar_number].NOAA_B = count IF count NE 0 THEN output[ar_number].NOAA_B_flares = fluxes(b_flares) c_flares = where (classes eq 'C',count) output[ar_number].NOAA_C = count IF count NE 0 THEN output[ar_number].NOAA_C_flares = fluxes(c_flares) m_flares = where (classes eq 'M',count) output[ar_number].NOAA_M = count IF count NE 0 THEN output[ar_number].NOAA_M_flares = fluxes(m_flares) x_flares = where (classes eq 'X',count) output[ar_number].NOAA_X = count IF count NE 0 THEN output[ar_number].NOAA_X_flares = fluxes(x_flares) ENDIF no_flares: ;gather other AR information ;this is messy as it involves a 2nd call to get_nar ;(already called inside flare_stats) ;need to clean this up IF noaa_reg NE -1 THEN BEGIN nar = get_nar2( map.time,/quiet,err=err,/unique ) IF err NE '' THEN BEGIN output[ar_number].Message=[output[ar_number].Message + 'Cannot get NAR info '] GOTO, cannot_get_NAR_info ENDIF nar = drot_nar2( nar, map.time, err=err) IF err NE '' THEN BEGIN output[ar_number].Message=[output[ar_number].Message + 'Cannot get NAR info '] GOTO, cannot_get_NAR_info ENDIF output[ar_number].NOAA_Size = nar(noaa_identified).area output[ar_number].NOAA_Class = string(nar(noaa_identified).ST$MAG_TYPE) output[ar_number].NOAA_Mac = string(nar(noaa_identified).ST$MACINTOSH) output[ar_number].NOAA_Posn_HG = nar(noaa_identified).location output[ar_number].NOAA_Posn_HC = [nar(noaa_identified).X, nar(noaa_identified).Y] output[ar_number].NOAA_Posn_EACP = xypos(*, noaa_identified) ENDIF cannot_get_NAR_info: new_filled=filled new_filled(where(filled LT -5.))=-10*(ar_number+1) ;so now new_filled is the AR pixel indices. ;*************************************************** ;these lines are a throwback to the original routine ;where I wanted to look at the contours and regions ;the contour map is no longer used here ;but the regmap is ;not do overwrite pixels which already belong to a previous AR contmap_notAR = where(contmap(0>(x_p-win_size):(sz_x-1)<(x_p+win_size-1),$ 0>(y_p-win_size):(sz_y-1)<(y_p+win_size-1)) EQ 0) ;move back into full-map array size subcont = contmap(0>(x_p-win_size):(sz_x-1)<(x_p+win_size-1),$ 0>(y_p-win_size):(sz_y-1)<(y_p+win_size-1)) subcont[ contmap_notAR ] = cor_cont[contmap_notAR] contmap(0>(x_p-win_size):(sz_x-1)<(x_p+win_size-1),$ 0>(y_p-win_size):(sz_y-1)<(y_p+win_size-1)) = subcont ;not do overwrite pixels which already belong to a previous AR regmap_notAR = where(regmap(0>(x_p-win_size):(sz_x-1)<(x_p+win_size-1),$ 0>(y_p-win_size):(sz_y-1)<(y_p+win_size-1)) GT -10) ;move back into full-map array size subreg = regmap(0>(x_p-win_size):(sz_x-1)<(x_p+win_size-1),$ 0>(y_p-win_size):(sz_y-1)<(y_p+win_size-1)) subreg[ regmap_notAR ] = new_filled[regmap_notAR] regmap(0>(x_p-win_size):(sz_x-1)<(x_p+win_size-1),$ 0>(y_p-win_size):(sz_y-1)<(y_p+win_size-1)) = subreg ;******************************************************** ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;define the data window for the output ;then flag all pixels which have already been ;accounted for in previous active regions subreg2 = copymap(0>(x_p-win_size):(sz_x-1)<(x_p+win_size-1),$ 0>(y_p-win_size):(sz_y-1)<(y_p+win_size-1)) prev_ar=fltarr(2.*win_size,2.*win_size) nonpixels = where(submap eq 0,count_np) IF count_np NE 0 THEN prev_ar(nonpixels) = 1 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;identify which pixels belong to this region filled_pixels=where(regmap eq -10*(ar_number+1)) ;identify positive pixels pos_pixels_num = where(carmap(filled_pixels) gt 0, count_pos) IF count_pos NE 0 THEN pos_pixels=filled_pixels(pos_pixels_num) ELSE pos_pixels =0 ;identify negative pixels neg_pixels_num = where(carmap(filled_pixels) lt 0, count_neg) IF count_neg NE 0 THEN neg_pixels=filled_pixels(neg_pixels_num) ELSE neg_pixels =0 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;now calculate stats for this AR ar_pos_size=N_ELEMENTS(pos_pixels)*1.4*1.4;conversion factor to get this into Mm^2 ar_neg_size=N_ELEMENTS(neg_pixels)*1.4*1.4;conversion factor to get this into Mm^2 ar_size=ar_pos_size+ar_neg_size ar_pos_fx=TOTAL(carmap(pos_pixels))*ar_pos_size/(10.^5.);*10^(15) to get into Mx ar_neg_fx=TOTAL(carmap(neg_pixels))*ar_neg_size/(10.^5.);*10^(15) to get into Mx ar_fx=ar_pos_fx+ar_neg_fx ar_fx_abs=ar_pos_fx+ABS(ar_neg_fx) contour_gaps=where(filled_contour eq 0) mix=fltarr(N_ELEMENTS(contour_gaps)-1) for i=0,N_ELEMENTS(mix)-1 do begin & $ start_contour=contour_gaps(i)+1 & $ end_contour=contour_gaps(i+1)-1 & $ if end_contour LT start_contour THEN GOTO, contour_skip if min(submap(filled_contour(start_contour:end_contour))) LT 0 AND $ max(submap(filled_contour(start_contour:end_contour))) GT 0 THEN $ mix(i)=1 & $ contour_skip: endfor num_mixed=N_ELEMENTS(where(mix eq 1)) IF count_pos eq 0 THEN BEGIN ar_pos_size = 0.0 ar_pos_fx = 0.0 ar_fx= ar_neg_fx ar_fx_abs = ABS(ar_neg_fx) mix(*) = 0.0 num_mixed=0.0 ENDIF IF count_neg eq 0 THEN BEGIN ar_neg_size = 0.0 ar_neg_fx= 0.0 ar_fx= ar_pos_fx ar_fx_abs = ABS(ar_pos_fx) mix(*) = 0.0 num_mixed=0.0 ENDIF output[ar_number].ARSE_data[ 0:(size(submap))(1)-1, 0: (size(submap))(2)-1 ] = subreg2 output[ar_number].ARSE_contour[ 0:(size(filled))(1)-1, 0:(size(filled))(2)-1 ] = filled output[ar_number].ARSE_nonpixels = prev_ar output[ar_number].ARSE_num = FIX(ar_number) output[ar_number].ARSE_posn_EACP = [round(x_p(0)),round(y_p(0))] output[ar_number].ARSE_Pos_Size = round(ar_pos_size) output[ar_number].ARSE_Neg_Size = round(ar_neg_size) output[ar_number].ARSE_Total_Size = round(ar_size) output[ar_number].ARSE_Pos_Flux = round(ar_pos_fx) output[ar_number].ARSE_Neg_Flux = round(ar_neg_fx) output[ar_number].ARSE_Total_Flux = round(ar_fx) output[ar_number].ARSE_Abs_Flux = round(ar_fx_abs) output[ar_number].ARSE_Mixing = [round(num_mixed),N_ELEMENTS(mix)] ;finished doing stats ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;need to remake carmap with this AR removed carmap(where(regmap NE 0))=0 ;now increase the ar_number and repeat the process ar_number=ar_number+1. IF (ar_number LT 49) THEN GOTO, find_next_region no_regions_remaining: ;********************** ;different selection criteria ;(1) once max(field) get below threshold, (-) arbitrary ;max_field=max(abs(carmap)) ;IF max_field GT 500 THEN GOTO, find_next_region ; ;(2) DEFAULT: ;if image cannot be contoured, stop. (+) not arbitrary (-) trouble with cosmic hits ;need to move 'GOTO, find_next_region' to before cannot_contour_image: and no_pixels_in_contour: ;also need to add in ;output=output(0:ar_number-1) at the bottom ; ;(3) once N active regions have been selected, (-) arbitrary ;IF (ar_number LT 15) THEN GOTO, find_next_region ;can act as a 'ridiculous outoput' switch - e.g., more than 50 regions is plain stupid ;*********************** IF ar_number eq 0 AND max (output[ar_number].ARSE_data) eq 0 THEN BEGIN err = 'No active regions selected' RETURN, -1 ENDIF IF ar_number eq 0 THEN RETURN, output(0) ;STOP RETURN, output(0:ar_number-1) END