code: def _numdiff_2d(func, x, y, dx=0, dy=0, eps=1e-8):
      if dx == 0 and dy == 0:
          return func(x, y)
      elif dx == 1 and dy == 0:
          return (func(x + eps, y) - func(x - eps, y)) / (2 * eps)
      elif dx == 0 and dy == 1:
          return (func(x, y + eps) - func(x, y - eps)) / (2 * eps)
      elif dx == 1 and dy == 1:
          return (
              func(x + eps, y + eps)
              - func(x - eps, y + eps)
              - func(x + eps, y - eps)
              + func(x - eps, y - eps)
          ) / (2 * eps) ** 2
      else:
          raise ValueError("invalid derivative order")
  
imports: import numpy
  from numpy.testing import assert_equal, assert_almost_equal, assert_array_equal, assert_array_almost_equal, assert_allclose, suppress_warnings
  from pytest import raises
  from numpy import array, diff, linspace, meshgrid, ones, pi, shape
  from scipy.interpolate.fitpack import bisplrep, bisplev
  from scipy.interpolate.fitpack2 import UnivariateSpline, LSQUnivariateSpline, InterpolatedUnivariateSpline, LSQBivariateSpline, SmoothBivariateSpline, RectBivariateSpline, LSQSphereBivariateSpline, SmoothSphereBivariateSpline, RectSphereBivariateSpline
label_desc: Use image node_id 1 for calling a global function with example usage: _numdiff_2d(func, x, y, dx, dy, eps) and returns: func, unknown, unknown, unknown