code: def sky2px(wcs,ra,dec,dra,ddec,cell, beam):
      
      dra = beam if dra<beam else dra # assume every source is at least as large as the psf
      ddec = beam if ddec<beam else ddec
      offsetDec = int((ddec/2.)/cell)
      offsetRA = int((dra/2.)/cell)
      if offsetDec%2==1:
          offsetDec += 1
      if offsetRA%2==1:
          offsetRA += 1
  
      raPix,decPix = map(int, wcs.wcs2pix(ra,dec))
      return np.array([raPix-offsetRA,raPix+offsetRA,decPix-offsetDec,decPix+offsetDec])
docs: convert a sky region to pixel positions