We carry out a nonlinear joint inversion of P wave first-motion polarities and coseismic surface displacement data of the 1908 Messina earthquake. We model the earthquake using a single planar fault: Slip is at first assumed to be uniform across the whole fault, then independent in a small set of coplanar subfaults, and finally smoothly variable across the fault. The first two steps are accomplished using a global minimization technique. The main features of the retrieved model are very robust and independent of the seismic velocity profile, of the presence of questionable bench marks, and of the assumed residual distribution. The along-strike component of slip is about half the along-dip component, in agreement with the direction of the extensional stress axis retrieved from geological observations. Surface vertical displacement is consistent with tsunami data and with the morphology of the Messina Straits.