;+ ; $Id: crop_by_coord.pro,v 1.10 2013/01/24 13:51:21 mcnutt Exp $ ; Project : STEREO - SECCHI ; ; Name : crop_by_coord ; ; Purpose : Extract a subfield from solar disk input image using the WCS. ; ; Explanation : Uses input coordinates to determine location of subfield to return. ; If location is behind limb AND /NOINTERP, center on point on limb nearest ; to input. ; ; Use : subfield = crop_by_coord(image, hdr, LLcorner, URcorner, outsize ) ; ; Inputs : image 2-d array ; hdr fits header or WCS header for image ; LLcorner 2-elem array with Carrington coordinate of of lower left corner of subfield ; -or- scalar CRLN of center of subfield ; URcorner 2-elem array with Carrington coordinate of of upper right corner of subfield ; -or- scalar CRLT of center of subfield ; (2nd value of URcorner is ignored; actual value determined by outsize) ; imsizeout scalar or 2-elem array specifying size of output array ; ; Optional Inputs : ; ; Outputs : FLTARR [outsize,outsize] or [outsize[0],outsize[1]] ; hdr fits header or WCS header for image (values updated) ; ; Keywords : ; /NOINTERP Do not project onto regular grid ; REF_HDR= IF /NOINTERP, use this to match platescale ; Calls : ; ; Category : Image Processing Display coordinates subfield wcs ; ; Prev. Hist. : None. ; ; Written : Nathan Rich, 2011 ; $Log: crop_by_coord.pro,v $ ; Revision 1.10 2013/01/24 13:51:21 mcnutt ; handles longitude crossing 360 ; ; Revision 1.9 2012/08/23 12:35:22 mcnutt ; set hdr crpixs vals based on crop if nointerp ; ; Revision 1.8 2012/08/22 19:52:32 nathan ; fix clt bug ; ; Revision 1.7 2012/08/21 22:02:50 nathan ; If location is behind limb AND /NOINTERP, center on point on limb nearest ; to input. ; ; Revision 1.6 2011/12/21 19:57:24 nathan ; correct x and y in computing lat and lng of default output ; ; Revision 1.5 2011/12/05 15:45:29 nathan ; fix comments ; ; Revision 1.4 2011/12/02 23:31:31 nathan ; allow subfi to go off edge of input; cancel change input hdr ; ; Revision 1.3 2011/12/02 22:55:16 nathan ; print info; return changed header (/nointerp only) ; ; Revision 1.2 2011/12/02 22:39:49 nathan ; allow input of center for either method ; ; Revision 1.1 2011/12/02 17:43:43 nathan ; works on individual images, header unchanged ; function crop_by_coord, img0, hdr0, llc, urc0, imsizeout, NOINTERP=nointerp, _EXTRA=_extra ; ;function euvi_heliographic, header, pixels, carrington=carrington, WCS=wcs ; ; Generate a WCS structure from the header, and extract the coordinates. Save ; the dimensions for later. ; IF tag_exist(hdr0,'WCSNAME') THEN wcs=hdr0 ELSE wcs = fitshead2wcs(hdr0) szp=size(img0) IF n_elements(imsizeout) LT 2 THEN imsizeout=[imsizeout,imsizeout] subim=fltarr(imsizeout[0],imsizeout[1]) urc=urc0 IF keyword_set(NOINTERP) THEN BEGIN ; In this case, take a subfield of correct size with no interpolation. ; ; ; IF urc[0] lt llc[0] THEN c360=360 ELSE c360=0 IF n_elements(llc) EQ 1 THEN cln0=llc ELSE cln0=llc[0]+ (urc[0]-llc[0]+c360)/2. IF cln0 gt 360 then cln0=cln0-360 IF n_elements(urc) EQ 1 THEN clt0=urc ELSE clt0=llc[1]+ (urc[1]-llc[1])/2. hdr=hdr0 img=scc_normalize(img0,hdr,interp=2,/norotate,_EXTRA=_extra) wcs=fitshead2wcs(hdr) WCS_CONVERT_to_COORD, wcs, hpcc, 'HG', cln0, clt0, /CARRINGTON clt=clt0[0] ; Check if cln,clt is visible. wf=where(finite(hpcc),nwf) IF nwf LT 2 THEN BEGIN print,'point is behind limb; centering on closest point on limb.' ; figure out which side point is in clnobs=hdr.crln_obs clnl=clnobs-90 clnr=clnobs+90 IF clnl LT 0 THEN clnl=clnl+360 IF clnr GT 360 THEN clnr=clnr-360 IF (cln0-clnr LT 90 and cln0-clnr GT 0) or (cln0-clnr LT -270) THEN BEGIN cln=clnr REPEAT BEGIN WCS_CONVERT_to_COORD, wcs, hpcc, 'HG', cln, clt, /CARRINGTON cln=cln-1 wf=where(finite(hpcc),nwf) ENDREP UNTIL nwf EQ 2 cln=cln+1 ENDIF ELSE BEGIN cln=clnl REPEAT BEGIN WCS_CONVERT_to_COORD, wcs, hpcc, 'HG', cln, clt, /CARRINGTON cln=cln+1 wf=where(finite(hpcc),nwf) ENDREP UNTIL nwf EQ 2 cln=cln-1 ENDELSE ENDIF ELSE BEGIN cln=cln0 ENDELSE IF ~(finite(cln)) THEN BEGIN message,'Chosen coordinate not valid for image; returning zeros.',/info wait,2 return,subim ENDIF cpix=wcs_get_pixel(wcs,hpcc) print,'Extracting ',trim(imsizeout[0]*hdr.cdelt1),'x',trim(imsizeout[1]*hdr.cdelt1),' arcsec box centered on [CRLN,CRLT]: ',trim(round(cln)),',',trim(round(clt)) ; If any coords are off the edge of the image, they should be zero in output. x1=cpix[0]-imsizeout[0]/2 x2=cpix[0]+imsizeout[0]/2-1 y1=cpix[1]-imsizeout[1]/2 y2=cpix[1]+imsizeout[1]/2-1 sx1=(0-x1)>0 sy1=(0-y1)>0 x1=x1>0 y1=y1>0 x2=x2<(hdr.naxis1-1) y2=y2<(hdr.naxis2-1) subim[sx1,sy1]=img[x1:x2, y1:y2] hdr.crpix1=hdr.crpix1-x1 hdr.crpix2=hdr.crpix2-y1 hdr.xcen=hdr.crpix1 hdr.ycen=hdr.crpix2 hdr.crpix1A=hdr.crpix1 hdr.crpix2A=hdr.crpix2 hdr0=hdr ENDIF ELSE BEGIN IF n_elements(llc) EQ 1 THEN BEGIN ; compute LL and UR corners ;RSUN=(6.95508e8 * 648d3 / !dpi ) / wcs.position.dsun_obs ;pscale1=60*wcs.cdelt[0]/rsun ; factor / solar radius in pix , hg arcsec -> carr deg pscale1=0.15 ; approx for 1024 llln=llc-pscale1*imsizeout[0]/2 lllt=urc-pscale1*imsizeout[1]/2 urln=llc+pscale1*(imsizeout[0]-1)/2 llc=[llln,lllt] urc=[urln,0.] ENDIF ; Figure out Lat of UR corner IF urc[0] lt llc[0] THEN c360=360 ELSE c360=0 pscale=(urc[0]-llc[0]+c360)/float(imsizeout[0]-1) urc[1]=llc[1]+pscale*(imsizeout[1]-1) print,'Extracting [CRLN1,CRLN2,CRLT1,CRLT2] ',string(llc[0],urc[0],llc[1],urc[1],format='(f8.3,f8.3,f8.3,f8.3)') lat=lgen(imsizeout[1],[llc[1],urc[1]])##replicate(1,imsizeout[0]) lng=replicate(1,imsizeout[1])##lgen(imsizeout[0],[llc[0],urc[0]+c360]) if c360 gt 0 then lng(where(lng gt 360))=lng((where(lng gt 360)))-360 ; Convert to HPC coordinates WCS_CONVERT_to_COORD, wcs, hpccoord, 'HG', lng, lat, /CARRINGTON pixels=wcs_get_pixel(wcs,hpccoord) print,'Interpolating image...' subim=reform(interpolate(float(img0),pixels[0,*,*],pixels[1,*,*],cubic=-0.5,missing=0.)) ENDELSE return,subim end