In practical applications of area-to-point spatial interpolation, inequality constraints, such as non-negativity, or more general constraints on the maximum and/or minimum allowable value of the resulting predictions, should be taken into account. The geostatistical framework proposed in this paper deals with area-to-point interpolation problems under such constraints, while: (i) explicitly accounting for support differences between sample data and unknown values, (ii) guaranteeing coherent predictions, and (iii) providing a measure of reliability for the resulting predictions. The analogy between the dual form of area-to-point interpolation and a spline allows to solve constrained area-to-point interpolation problems via a constrained quadratic minimization algorithm, after accounting for the following three issues: (i) equality and inequality constraints could be applied to different supports, and such support differences should be considered explicitly in the problem formulation, (ii) if inequality constraints are enforced on the entire set of points discretizing the areal data, it is impossible to obtain a solution of the quadratic programming problem, and (iii) the uniqueness and existence of the solution has to be diagnosed. In this work, stable and efficient computation of point predictions is achieved through the following two steps: (i) initial prediction at all locations via unconstrained area-to-point interpolation, and (ii) constrained area-to-point interpolation with inequality information only at those points whose initial predicted values violate the inequality constraints. Last, the application of the proposed method to area-to-point spatial interpolation with inequality constraints in one and two dimensions is demonstrated using realistically simulated data.