Gaussian processes (GPs) are instrumental in modeling spatial processes, offering precise interpolation and prediction capabilities across fields such as environmental science and biology. Recently, there has been growing interest in extending GPs to infer spatial derivatives, which are vital for analyzing spatial dynamics and detecting subtle changes in data patterns. Despite their utility, traditional GPs suffer from computational inefficiencies, due to the cubic scaling with the number of spatial locations. Fortunately, the computational challenge has spurred extensive research on scalable GP methods. However, these scalable approaches do not directly accommodate the inference of derivative processes. A straightforward approach is to use scalable GP models followed by finite-difference methods, known as the plug-in estimator. This approach, while intuitive, suffers from sensitivity to parameter choices, and the approximate gradient may not be a valid GP, leading to compromised inference. To bridge this gap, we introduce the Nearest-Neighbor Derivative Process (NNDP), an innovative framework that models the spatial processes and their derivatives within a single scalable GP model. NNDP significantly reduces the computational time complexity from $O(n^3)$ to $O(n)$, making it feasible for large datasets. We provide various theoretical supports for NNDP and demonstrate its effectiveness through extensive simulations and real data analysis.