Abstract: The Wilson depression is the difference in geometric height of the layer of unit continuum optical depth between the sunspot umbra and the quiet Sun. Measuring the Wilson depression is important for understanding the geometry of sunspots. Current methods suffer from systematic effects or need to make assumptions on the geometry of the magnetic field. This leads to large systematic uncertainties of the derived Wilson depressions. Here we present a method for deriving the Wilson depression that only requires the information about the magnetic field that are accessible by spectropolarimetry and that does not rely on assumptions on the geometry of sunspots or on its magnetic field. Our method is based on minimizing the divergence of the magnetic field vector derived from spectropolarimetric observations. We focus on large spatial scales only in order to reduce the number of free parameters. We test the performance of our method using synthetic Hinode data derived from two sunspot simulations. We find that the maximum and the umbral averaged Wilson depression for both spots determined with our method typically lies within 100 km of the true value obtained from the simulations. In addition, we apply the method to spots from the Hinode sunspot database at MPS. The derived Wilson depressions (500-700 km) are consistent with results typically obtained from the Wilson effect. In our sample, larger spots with a stronger magnetic field exhibit a higher Wilson depression than smaller spots.