title: How to obtain 2D Cutout of an image from a SkyCoord position?
body: <p>I am following the example from Astropy docs for 2D Cutout. </p>
  
  <p>The header of my FITS file:</p>
  
  <pre><code>SIMPLE  =                    T / file does conform to FITS standard             
  BITPIX  =                  -32 / number of bits per data pixel                  
  NAXIS   =                    3 / number of data axes                            
  NAXIS1  =                  512 / length of data axis 1                          
  NAXIS2  =                  512 / length of data axis 2                          
  NAXIS3  =                    3 / length of data axis 3                          
  EXTEND  =                    T / FITS dataset may contain extensions            
  COMMENT   FITS (Flexible Image Transport System) format is defined in 'Astronomy
  COMMENT   and Astrophysics', volume 376, page 359; bibcode: 2001A&amp;A...376..359H 
  SURVEY  = 'DECaLS  '                                                            
  VERSION = 'DR8-south'                                                           
  IMAGETYP= 'image   '                                                            
  BANDS   = 'grz     '                                                            
  BAND0   = 'g       '                                                            
  BAND1   = 'r       '                                                            
  BAND2   = 'z       '                                                            
  CTYPE1  = 'RA---TAN'           / TANgent plane                                  
  CTYPE2  = 'DEC--TAN'           / TANgent plane                                  
  CRVAL1  =            186.11382 / Reference RA                                   
  CRVAL2  =           0.15285422 / Reference Dec                                  
  CRPIX1  =                256.5 / Reference x                                    
  CRPIX2  =                256.5 / Reference y                                    
  CD1_1   = -7.27777777777778E-05 / CD matrix                                     
  CD1_2   =                   0. / CD matrix                                      
  CD2_1   =                   0. / CD matrix                                      
  CD2_2   = 7.27777777777778E-05 / CD matrix                                      
  IMAGEW  =                 512. / Image width                                    
  IMAGEH  =                 512. / Image height                                   
  </code></pre>
  
  <p>So far what I have tried :</p>
  
  <pre><code>from astropy.coordinates import SkyCoord
  from astropy.wcs import WCS
  
  position = SkyCoord(hdu[0].header['CRVAL1']*u.deg,hdu[0].header['CRVAL2']*u.deg)
  size = 200*u.pixel
  
  wcs1 = WCS(hdu[0].header)
  
  cutout = Cutout2D(hdu[0].data[0], position ,size, wcs = wcs1 )
  </code></pre>
  
  <p>I run into error in the last line. 
  Error : </p>
  
  <pre><code>---------------------------------------------------------------------------
  ValueError                                Traceback (most recent call last)
  &lt;ipython-input-142-7cc21a13e941&gt; in &lt;module&gt;
  ----&gt; 1 cutout = Cutout2D(hdu[0].data[0], position ,size, wcs = wcs1 )
  
  /Applications/anaconda3/lib/python3.7/site-packages/astropy/nddata/utils.py in __init__(self, data, position, size, wcs, mode, fill_value, copy)
      714         if wcs is not None:
      715             self.wcs = deepcopy(wcs)
  --&gt; 716             self.wcs.wcs.crpix -= self._origin_original_true
      717             self.wcs.array_shape = self.data.shape
      718             if wcs.sip is not None:
  
  ValueError: operands could not be broadcast together with shapes (3,) (2,) (3,) 
  </code></pre>
  
  <p>My guess is it is because naxis =3 in my file and the documentation assumes naxis = 2. Though I am not sure if that is the actual issue here. Can anybody help fix this error?</p>