Is your original function only on an (x',y') grid (e.g. from a finite-difference calculation)? If so, you are probably stuck. (You could interpolate as @Oscar_Smith suggests, but you generally can’t get high-order accuracy interpolating from a regular grid … unless it is periodic.) You could use Romberg integration on an equispaced grid, however.