We present two methods based on functional principal component analysis (FPCA) for the estimation of smooth derivatives of a sample of random functions, which are observed in a more than one-dimensional domain. We apply eigenvalue decomposition to a) the dual covariance matrix of the derivatives, and b) the dual covariance matrix of the observed curves. To handle noisy data from discrete observations, we rely on local polynomial regressions. If curves are contained in a finite-dimensional function space, the second method performs better asymptotically. We apply our methodology in a simulation and empirical study, in which we estimate state price density (SPD) surfaces from call option prices. We identify three main components, which can be interpreted as volatility, skewness and tail factors. We also find evidence for term structure variation.