Radial basis functions have gained popularity for many applications including numerical solution of partial differential equations, image processing, and machine learning. For these applications it is useful to have an algorithm which detects edges or sharp gradients and is based on the underlying basis functions. In our previous research, we proposed an iterative adaptive multiquadric radial basis function method for the detection of local jump discontinuities in one-dimensional problems. The iterative edge detection method is based on the observation that the absolute values of the expansion coefficients of multiquadric radial basis function approximation grow exponentially in the presence of a local jump discontinuity with fixed shape parameters but grow only linearly with vanishing shape parameters. The different growth rate allows us to accurately detect edges in the radial basis function approximation. In this work, we extend the one-dimensional iterative edge detection method to two-dimensional problems. We consider two approaches: the dimension-by-dimension technique and the global extension approach. In both cases, we use a rescaling method to avoid ill-conditioning of the interpolation matrix. The global extension approach is less efficient than the dimension-by-dimension approach, but is applicable to truly scattered two-dimensional points, whereas the dimension-by-dimension approach requires tensor product grids. Numerical examples using both approaches demonstrate that the two-dimensional iterative adaptive radial basis function method yields accurate results.