The gully boundary, which distinguishes gully from non‐gully areas, is a significant geomorphologic feature for research on gully development and gully erosion. This study presents a new method based on multidirectional hill‐shading maps, which identify the gully (i.e., shadow area) from inter‐gully (i.e., non‐shadow area). These shadows obtained from various illumination azimuths are merged; consequently, the border of the shadows, which is the gully boundary line, can be achieved. In this process, two key parameters, namely, altitude and azimuth of light, affect the accuracy of gully boundary extraction. The experiments in Yaojiawan area of China show that the method of average median slope of all sampling profiles across the gully boundary is effective and practical for light altitude selection. Moreover, the six azimuths are sufficient for gully boundary extraction in a loess hilly area. The application in the Madigou area indicates the replicability and rationality of this method. A comparison with the positive‐negative terrain and slope variation method confirms a higher accuracy of gully boundary extraction by the proposed method in terms of visual interpretation, length, and contour‐matching difference with reference to manually digitalized results. Accuracy assessment indicates that the proposed method is applicable for gully boundary extraction based on high‐resolution DEMs.