lorentzian - python fit power law

Sé que scipy curve_fit puede mejorar

Question

Estoy usando python / numpy / scipy para implementar este algoritmo para alinear dos modelos digitales de elevación (DEM) en función del aspecto del terreno y la pendiente:

"Co-registro y correcciones de sesgo de los conjuntos de datos de elevación de satélites para cuantificar el cambio de espesor del glaciar", C. Nuth y A. Kääb, doi: 10.5194 / tc-5-271-2011

Tengo todo un marco configurado, pero la calidad del ajuste proporcionado por scipy.optimize.curve_fit es deficiente.

def f(x, a, b, c):
    y = a * numpy.cos(numpy.deg2rad(b-x)) + c
    return y

def compute_offset(dh, slope, aspect):
    import scipy.optimize as optimization

    idx = random.sample(range(dh.compressed().size), 10000)
    xdata = numpy.array(aspect.compressed()[idx], float)
    ydata = numpy.array((dh/numpy.tan(numpy.deg2rad(slope))).compressed()[idx], float)

    #Generate synthetic data to test curve_fit
    #xdata = numpy.arange(0,360,0.01)
    #ydata = f(xdata, 20.0, 130.0, -3.0) + 20*numpy.random.normal(size=len(xdata))

    print xdata
    print ydata

    x0 = numpy.array([0.0, 0.0, 0.0])

    fit = optimization.curve_fit(f, xdata, ydata, x0)[0]
    #optimization.leastsq(f, x0[:], args=(xdata, ydata))
    genplot(xdata, ydata, fit)
    return fit

def genplot(x, y, fit):
    a = (numpy.arange(0,360))
    f_a = f(a, fit[0], fit[1], fit[2])
    idx = random.sample(range(x.size), 10000)
    plt.figure()
    plt.xlabel(''Aspect (deg)'')
    plt.ylabel(''dh/tan(slope) (m)'')
    plt.plot(x[idx], y[idx], ''r.'')
    plt.axhline(color=''k'')
    plt.plot(a, f_a, ''b'')
    plt.ylim(-80,80)
    plt.show()

#Input DEMs
dem1_fn = sys.argv[1]
dem2_fn = sys.argv[2]

dem1_ds = gdal.Open(dem1_fn, gdal.GA_ReadOnly)
dem2_ds = gdal.Open(dem2_fn, gdal.GA_ReadOnly)

#Extract band 1 from each dataset as masked array using internal nodata value
dem1 = getperc_new.gdal_getma(dem1_ds, 1)
dem2 = getperc_new.gdal_getma(dem2_ds, 1)

#Produce slope and aspect maps using gdaldem and load into masked arrays
dem1_slope = gdaldem_slope(dem1_fn)
dem1_aspect = gdaldem_aspect(dem1_fn)

#Compute common mask and apply to all products
common_mask = dem1.mask + dem2.mask + dem1_slope.mask + dem1_aspect.mask

diff_euler = numpy.ma.array(dem2-dem1, mask=common_mask)
dem1_slope.__setmask__(common_mask)
dem1_aspect.__setmask__(common_mask)

#Compute relationship between elevation difference, slope and aspect
fit = compute_offset(diff_euler, dem1_slope, dem1_aspect)
print fit

Aquí está el ajuste para mis datos, que inicialmente consisten en ~ 2 millones de puntos, pero he tomado muestras al azar para fines de prueba / trazado:

[-14.9639559 216.01093596 -41.96806735]

Hay muchos datos para un buen ajuste, pero el resultado de curve_fit es pobre. Cuando corro con datos sintéticos, obtengo un buen ajuste:

parámetros de entrada originales [20.0, 130.0, -3.0]

resultado de curve_fit [-19.66719631 -49.6673076 -3.12198723]

No estoy seguro si esto tiene algo que ver con el uso de matrices enmascaradas, una limitación de curve_fit, o si simplemente estoy pasando por alto algo simple. Gracias por cualquier sugerencia.

=======================

Editar 9/4/13 16:30 PDT

Según lo sugerido por @Evert y otros, el problema definitivamente estaba relacionado con valores atípicos. Pude obtener un ajuste mucho mejor después de eliminar los valores atípicos. Mirando mi código anterior, parece que calculé la desviación absoluta media para cada aspecto, luego eliminé todo lo que estuviera fuera de 2 * enojado antes de ajustarlo.

Genere algunas parcelas adicionales en noviembre de 2012:

Pero volviendo a ver esto, estoy casi seguro de que se generaron para diferentes datos de entrada. Es todo lo que puedo encontrar en este momento, así que los incluyo aquí como un ejemplo de un caso de muestreo sesgado. Este método para la alineación DEM está destinado a fallar en casos como estos, y no tiene nada que ver con las capacidades de ajuste de curvas de scipy.

Terminé desarrollando un enfoque diferente para la alineación que implicaba la correlación cruzada normalizada, el refinamiento subpíxel y la eliminación del desplazamiento vertical para dos matrices numpy 2D enmascaradas. Es más rápido y proporciona consistentemente mejores resultados. Aunque incluso ese enfoque ha sido reemplazado por una herramienta Iterative Closest Point (ICP) (pc_align) desarrollada por Oleg Alexandrov como parte de la Ames Stereo Pipeline de la NASA .

Gracias por todas sus respuestas y me disculpo por abandonar esta pregunta.

Answer 1

Si solo intenta obtener una onda sinusoidal con compensación de fase, no necesita un ajuste no lineal.

Puedes reemplazar ese a * sin(x - b) + c por a * sin(x) + b * cos(x) + c , porque cualquier seno con un desplazamiento se puede escribir como una combinación apropiada de un seno y un coseno ("Adición de phaser", como en una transformada de Fourier).

Si eso da el mismo resultado, entonces ese problema no es el ajuste "no lineal".