function ds = cosmo_surface_dataset(fn, varargin)
cosmo structjoin hdr
cosmo surficial neighborhood hdr