Here I develop a cross-correlation approach to estimating heights of shoreline features, and apply the new method to paleo-shorelines of Pleistocene Lake Bonneville. I calculate 1st-derivative (slope) and 2nd-derivative (curvature) profiles from Digital Elevation Model (DEM) or Global Positioning System Real-Time Kinematic (GPS-RTK) measurements of elevation. I then cross-correlate pairs of profiles that have been shifted by various "lags," or shifts in elevation. The correlation coefficient (a normalized dot-product measure of similarity) is calculated as a function of lag within small (~40 m) windows centered at various elevations. The elevation and lag with the greatest correlation coefficient indicates the shoreline elevation at the reference profile and the change in shoreline height for the profile pair. I evaluate several different algorithms for deriving slope and curvature by examining closure of elevation lags across profile triples. I then model isostatic response to Lake Bonneville loading and unloading. I first model lakeshore uplift response to lake load removal assuming an elastic layer over an inviscid half-space. I obtain a best-fit comparison of predicted to observed shoreline heights for the Bonneville level with an elastic layer thickness, Te, of 25±2 km (at 95% confidence) when using only previously published shoreline elevation estimates. The best-fit for the Bonneville level when using these estimates plus 44 new estimates suggests a Te of 26±2 km. The best-fit model for the Provo level suggests Te of 17±3 km. For the Gilbert level, the response is insensitive to the assumed Te. I next model isostatic response to Bonneville loading and unloading assuming an elastic layer over a viscoelastic halfspace. This approach assumes constant parameters for the entire loading history, and yields a best-fit model with Te =70±5 km and viscosity ç=~2x1018 Pa s with 95% confidence ranging from ~1x1018 to ~5x1019 Pa s when only the previously published data are used. With the newer data added, the best-fit model has Te =58±2 km and ç ranging from ~1x1018 to ~1x1019 Pa s with 95% confidence. The 12-15 m weighted root-mean-square misfit to the best-fitting model is dominated by tectonic signals related to Basin-and-Range tectonics particularly seismic offsets of the Wasatch fault, and closely mimics the geological timescale pattern of basin-subsidence and range-uplift.