We observe $n$ inhomogeneous Poisson processes with covariates and aim at estimating their intensities. We assume that the intensity of each Poisson process is of the form $s (\cdot, x)$ where $x$ is the covariate and where $s$ is an unknown function. We propose a model selection approach where the models are used to approximate the multivariate function $s$. We show that our estimator satisfies an oracle-type inequality under very weak assumptions both on the intensities and the models. By using an Hellinger-type loss, we establish non-asymptotic risk bounds and specify them under several kind of assumptions on the target function $s$ such as being smooth or a product function. Besides, we show that our estimation procedure is robust with respect to these assumptions.