Nonparametric estimation of probability density functions, both marginal and joint densities, is a very useful tool in statistics. The kernel method is popular and applicable to dependent data, including time series and spatial data. But at least for the joint density, one has had to assume that data are observed at regular time intervals or on a regular grid in space. Though this is not very restrictive in the time series case, it often is in the spatial case. In fact, to a large degree it has precluded applications of nonparametric methods to spatial data because such data often are irregularly positioned over space. In this article, we propose nonparametric kernel estimators for both the marginal and in particular the joint probability density functions for nongridded spatial data. Large sample distributions of the proposed estimators are established under mild conditions, and a new framework of expanding-domain infill asymptotics is suggested to overcome the shortcomings of spatial asymptotics in the existing literature. A practical, reasonable selection of the bandwidths on the basis of cross-validation is also proposed. We demonstrate by both simulations and real data examples of moderate sample size that the proposed methodology is effective and useful in uncovering nonlinear spatial dependence for general, including non-Gaussian, distributions. Supplementary materials for this article are available online.