Let {(Yi, Xi), i∈ℤN} be a stationary real-valued (d+1)-dimensional spatial processes. Denote by x↦qp(x), p∈(0, 1), x∈ℝd, the spatial quantile regression function of order p, characterized by P{Yi≤qp(x)|Xi=x}=p. Assume that the process has been observed over an N-dimensional rectangular domain of the form $\mathcal{I}_{\mathbf{n}}:=\{\mathbf{i}=(i_{1},\ldots,i_{N})\in\mathbb{Z}^{N}\vert1\leq i_{k}\leq n_{k},k=1,\ldots,N\}$ , with n=(n1, …, nN)∈ℤN. We propose a local linear estimator of qp. That estimator extends to random fields with unspecified and possibly highly complex spatial dependence structure, the quantile regression methods considered in the context of independent samples or time series. Under mild regularity assumptions, we obtain a Bahadur representation for the estimators of qp and its first-order derivatives, from which we establish consistency and asymptotic normality. The spatial process is assumed to satisfy general mixing conditions, generalizing classical time series mixing concepts. The size of the rectangular domain $\mathcal{I}_{\mathbf{n}}$ is allowed to tend to infinity at different rates depending on the direction in ℤN (non-isotropic asymptotics). The method provides much richer information than the mean regression approach considered in most spatial modelling techniques.