def compute_offset_nuth(dh, slope, aspect, min_count=<NUM_LIT:100>, remove_outliers=True, plot=True):
import scipy.optimize as optimization<EOL>if dh.count() < min_count:<EOL><INDENT>sys.exit("<STR_LIT>")<EOL><DEDENT>if slope.count() < min_count:<EOL><INDENT>sys.exit("<STR_LIT>")<EOL><DEDENT>med_dh = malib.fast_median(dh)<EOL>med_slope = malib.fast_median(slope)<EOL>c_seed = (med_dh/np.tan(np.deg2rad(med_slope))) <EOL>x0 = np.array([<NUM_LIT:0.0>, <NUM_LIT:0.0>, c_seed])<EOL>print("<STR_LIT>")<EOL>common_mask = ~(malib.common_mask([dh, aspect, slope]))<EOL>xdata = aspect[common_mask].data<EOL>ydata = (dh[common_mask]/np.tan(np.deg2rad(slope[common_mask]))).data<EOL>print("<STR_LIT>")<EOL>print(ydata.size)<EOL>if remove_outliers:<EOL><INDENT>print("<STR_LIT>")<EOL>f = <NUM_LIT:3><EOL>sigma, u = (ydata.std(), ydata.mean())<EOL>rmin = u - f*sigma<EOL>rmax = u + f*sigma<EOL>print("<STR_LIT>" % (rmin, rmax))<EOL>idx = (ydata >= rmin) & (ydata <= rmax)<EOL>xdata = xdata[idx]<EOL>ydata = ydata[idx]<EOL>print(ydata.size)<EOL><DEDENT>nbins = <NUM_LIT><EOL>bin_range = (<NUM_LIT:0.>, <NUM_LIT>)<EOL>bin_width = <NUM_LIT:1.0><EOL>bin_count, bin_edges, bin_centers = malib.bin_stats(xdata, ydata, stat='<STR_LIT:count>', nbins=nbins, bin_range=bin_range)<EOL>bin_med, bin_edges, bin_centers = malib.bin_stats(xdata, ydata, stat='<STR_LIT>', nbins=nbins, bin_range=bin_range)<EOL>"""<STR_LIT>"""<EOL>min_bin_sample_count = <NUM_LIT:9><EOL>idx = (bin_count.filled(<NUM_LIT:0>) >= min_bin_sample_count) <EOL>bin_count = bin_count[idx].data<EOL>bin_med = bin_med[idx].data<EOL>bin_centers = bin_centers[idx]<EOL>fit = None<EOL>fit_fig = None<EOL>min_bin_count = <NUM_LIT> <EOL>bin_ptp = np.cos(np.radians(bin_centers)).ptp()<EOL>min_bin_ptp = <NUM_LIT:1.0> <EOL>if len(bin_med) >= min_bin_count and bin_ptp >= min_bin_ptp:<EOL><INDENT>print("<STR_LIT>")<EOL>fit = optimization.curve_fit(nuth_func, bin_centers, bin_med, x0)[<NUM_LIT:0>]<EOL>print(fit)<EOL>if plot:<EOL><INDENT>print("<STR_LIT>")<EOL>bin_idx = np.digitize(xdata, bin_edges)<EOL>output = []<EOL>for i in np.arange(<NUM_LIT:1>, len(bin_edges)):<EOL><INDENT>output.append(ydata[bin_idx==i])<EOL><DEDENT>lw = <NUM_LIT><EOL>whiskerprops={'<STR_LIT>':lw}<EOL>capprops={'<STR_LIT>':lw}<EOL>boxprops={'<STR_LIT>':'<STR_LIT:k>', '<STR_LIT>':<NUM_LIT:0>}<EOL>medianprops={'<STR_LIT>':'<STR_LIT:o>', '<STR_LIT>':<NUM_LIT:1>, '<STR_LIT>':'<STR_LIT:r>'}<EOL>fit_fig, ax = plt.subplots(figsize=(<NUM_LIT:6>,<NUM_LIT:6>))<EOL>widths = <NUM_LIT>*(bin_count/bin_count.max())<EOL>s=<NUM_LIT:3><EOL>bp = ax.boxplot(np.array(output)[idx][::s], positions=bin_centers[::s], widths=widths[::s], showfliers=False,patch_artist=True, boxprops=boxprops, whiskerprops=whiskerprops, capprops=capprops,medianprops=medianprops)<EOL>bin_ticks = [<NUM_LIT:0>, <NUM_LIT>, <NUM_LIT>, <NUM_LIT>, <NUM_LIT>, <NUM_LIT>, <NUM_LIT>, <NUM_LIT>, <NUM_LIT>]<EOL>ax.set_xticks(bin_ticks)<EOL>ax.set_xticklabels(bin_ticks)<EOL>"""<STR_LIT>"""<EOL>f_a = nuth_func(bin_centers, fit[<NUM_LIT:0>], fit[<NUM_LIT:1>], fit[<NUM_LIT:2>])<EOL>nuth_func_str = r'<STR_LIT>' % tuple(fit)<EOL>ax.plot(bin_centers, f_a, '<STR_LIT:b>', label=nuth_func_str)<EOL>ax.set_xlabel('<STR_LIT>')<EOL>ax.set_ylabel('<STR_LIT>')<EOL>ax.axhline(color='<STR_LIT>', linewidth=<NUM_LIT:0.5>)<EOL>ax.set_xlim(*bin_range)<EOL>ylim = ax.get_ylim()<EOL>abs_ylim = np.max(np.abs(ylim))<EOL>pad = <NUM_LIT:0><EOL>ylim = (-abs_ylim - pad, abs_ylim + pad)<EOL>minylim = (-<NUM_LIT:10>,<NUM_LIT:10>)<EOL>if ylim[<NUM_LIT:0>] > minylim[<NUM_LIT:0>]:<EOL><INDENT>ylim = minylim<EOL><DEDENT>ax.set_ylim(*ylim)<EOL>ax.legend(prop={'<STR_LIT:size>':<NUM_LIT:8>})<EOL><DEDENT><DEDENT>return fit, fit_fig<EOL>
Compute horizontal offset between input rasters using Nuth and Kaab [2011] (nuth) method