We propose a novel method for the analysis of spatially distributed data from an exponential family distribution, able to efficiently treat data occurring over irregularly shaped domains. We consider a generalized linear framework and extend the work of Sangalli et al. (2013) to distributions other than the Gaussian. In particular, we can handle all distributions within the exponential family, including binomial, Poisson and Gamma outcomes, hence leading to a very broad applicability of the proposed model. We maximize a penalized log-likelihood function. The roughness penalty term involves a suitable differential operator of the spatial field over the domain of interest. This maximization is done via a penalized iterative least square approach (see Wood (2006)). Covariate information can also be included in the model in a semi-parametric setting. The proposed models exploit advanced scientific computing techniques and specifically make use of the Finite Element Method, that provides a basis for piecewise polynomial surfaces and allows to impose boundary conditions on the space distribution of the probability. Finally, we extend theoretically the model to deal with data occurring on a two dimensional manifold.