We introduce a new variational estimator for the intensity function of an inhomogeneous spatial point process with points in the $d$-dimensional Euclidean space and observed within a bounded region. The variational estimator applies in a simple and general setting when the intensity function is assumed to be of log-linear form $β+{θ}^{\top}z(u)$ where $z$ is a spatial covariate function and the focus is on estimating ${θ}$. The variational estimator is very simple to implement and quicker than alternative estimation procedures. We establish its strong consistency and asymptotic normality. We also discuss its finite-sample properties in comparison with the maximum first order composite likelihood estimator when considering various inhomogeneous spatial point process models and dimensions as well as settings were $z$ is completely or only partially known.