Respuesta corta: Es un problema de precisión numérica: el xa , xb y xc Los coeficientes no se muestran en el archivo de salida 6S con la suficiente precisión para producir la respuesta correcta cuando se utiliza la fórmula. Es posible modificar el 6S para que imprima los coeficientes con una mayor precisión, lo que resuelve el problema.
Respuesta larga: Me intrigó esta pregunta, ya que he utilizado 6S durante mucho tiempo (y he escrito una interfaz de Python para ello) y nunca me había dado cuenta de esto... ¡extrañamente, nunca había pensado en comprobarlo!
En un principio pensé que podría tratarse de un problema con tu parametrización del 6S, pero probé algunas parametrizaciones diferentes y me encontré con el mismo problema: obtenía una reflectancia corregida atmosféricamente ligeramente diferente al pasar los coeficientes por la fórmula, en comparación con la reflectancia que arrojaba el modelo 6S directamente.
El manual del 6S es muy detallado, pero de alguna manera nunca parece responder a las preguntas que tengo - por ejemplo, no explica en ninguna parte cómo se calculan los tres coeficientes. Sin embargo, tiene un archivo de salida de ejemplo que incluye los resultados de la corrección atmosférica (ver la página final de Parte 1 del manual ). Esto incluye las siguientes salidas:
*******************************************************************************
* atmospheric correction result *
* ----------------------------- *
* input apparent reflectance : 0.100 *
* measured radiance [w/m2/sr/mic] : 38.529 *
* atmospherically corrected reflectance *
* Lambertian case : 0.22180 *
* BRDF case : 0.22180 *
* coefficients xa xb xc : 0.00685 0.03885 0.06835 *
* y=xa*(measured radiance)-xb; acr=y/(1.+xc*y) *
*******************************************************************************
Si se trabaja en el cálculo utilizando la fórmula dada se encuentra que el resultado del cálculo no coincide con la salida del 6S. Permítanme decirlo de nuevo: en el ejemplo proporcionado por los autores del 6S, ¡la salida del modelo y la fórmula no coinciden! No podía creer esto...
Así que me pregunté si la fórmula era algún tipo de ajuste de curva simple a algunas salidas del 6S, y por lo tanto se esperaría que tuviera un pequeño error en comparación con las salidas reales del modelo. Como ya he mencionado, el manual explica muchas cosas con gran detalle, pero no dice nada sobre el cálculo de estos coeficientes. Por suerte, el código fuente del 6S está disponible para descargar . Y lo que es menos conveniente, el código fuente está escrito en Fortran 77.
No soy en absoluto un experto en Fortran 77 (de hecho, nunca he escrito ningún código Fortran en la vida real), pero he indagado en el código para intentar averiguar cómo se calculan los coeficientes.
Si quieres seguirlo, el código para calcular los coeficientes comienza en la línea 3382 de main.f . Los coeficientes reales se fijan en las líneas 3393-3397:
xa=pi*sb/xmus/seb/tgasm/sutott/sdtott
xap=1./tgasm/sutott/sdtott
xb=ainr(1,1)/sutott/sdtott/tgasm
xb=ainr(1,1)/sutott/sdtott/tgasm
xc=sast
(extrañamente xb se fija dos veces, con el mismo valor, y otro coeficiente xap que parece que nunca se utiliza (¡no tengo ni idea de por qué!).
A partir de este código es bastante obvio que no se utiliza ningún algoritmo complicado de ajuste de curvas: los coeficientes son simplemente manipulaciones algebraicas de otras variables utilizadas en el modelo. Por ejemplo, xc se establece el valor de la variable sast que, con un poco de trabajo de detective, resulta ser el albedo esférico total (ver línea 3354). Puede comprobarlo en la salida del 6S: el valor de xc es siempre el mismo que el albedo esférico total que se muestra unas líneas más arriba en el archivo de salida. Del mismo modo, xb se calcula en función de diversas variables, entre ellas tgasm que es la transmitancia global total del gas y sdtott que es la transmitancia total de la dispersión hacia abajo, y así sucesivamente. (Estas variables son difíciles de descifrar, porque Fortran 77 tiene un límite de seis caracteres para los nombres de las variables, por lo que no son muy descriptivas).
Estaba perplejo en este punto, hasta que pensé en la precisión numérica. Me di cuenta de que el xa El coeficiente tiene un número de ceros después del punto decimal, y se pregunta si no hay suficientes cifras significativas para producir un resultado preciso al utilizar la fórmula. Resultó que este era el caso, pero voy a ir a través de cómo he alterado el código 6S para probar esto.
Línea 3439 de main.f se encarga de escribir los coeficientes en el archivo. Se compone de:
write(iwr, 944)xa,xb,xc
Esto le dice a Fortran que escriba la salida en el flujo de archivo/salida iwr utilizando el código de formato especificado en la línea 944 y escribir las tres variables xa , xb y xc . Si miramos la línea 944 (es decir, la línea a la que se le ha dado un número de línea Fortran de 944, que en realidad es la línea 3772 del archivo... ¡sólo para mantenerte alerta!
944 format(1h*,6x,40h coefficients xa xb xc :,
s 1x, 3(f8.5,1x),t79,1h*,/,1h*,6x,
s ' y=xa*(measured radiance)-xb; acr=y/(1.+xc*y)',
s t79,1h*,/,79(1h*))
Esta línea bastante complicada explica cómo formatear la salida. El bit clave es 3(f8.5,1x) que le dice a Fortran que escriba un número de punto flotante ( f ) con una anchura máxima de 8 caracteres y 5 decimales ( 8.5 ) seguido de un espacio ( 1x ), y repetirlo tres veces (el 3(...) ). Podemos modificar esto para imprimir más decimales - por ejemplo, yo lo cambié a 3(f10.8,1x) , lo que nos da 8 decimales. Si hacemos esto, entonces encontramos que la salida va hacia el * que están al final de cada línea, por lo que tenemos que modificar un poco el resto de la línea para reducir el número de espacios después del texto coefficients xa xb xc . La línea final, que funciona, tiene este aspecto:
944 format(1h*,6x,35h coefficients xa xb xc :,
s 1x, 3(f10.8,1x),t79,1h*,/,1h*,6x,
s ' y=xa*(measured radiance)-xb; acr=y/(1.+xc*y)',
s t79,1h*,/,79(1h*))
Si se modifica esta línea en main.f y recompilar el 6S, verá que su salida se parece a esto:
*******************************************************************************
* atmospheric correction result *
* ----------------------------- *
* input apparent reflectance : 0.485 *
* measured radiance [w/m2/sr/mic] : 240.000 *
* atmospherically corrected reflectance *
* Lambertian case : 0.45439 *
* BRDF case : 0.45439 *
* coefficients xa xb xc : 0.00297362 0.20291930 0.24282509 *
* y=xa*(measured radiance)-xb; acr=y/(1.+xc*y) *
*******************************************************************************
Si a continuación se aplica la fórmula, se comprobará que el resultado de la fórmula y el resultado del modelo coinciden, al menos hasta el número de decimales del resultado del modelo.
En mis pruebas de esto, obtuve lo siguiente para el código original 6S:
- Modelo: 0.4543900000
- Fórmula: 0.4537049078
- Perc Diff: 0.1507718536%
(la diferencia porcentual que yo obtenía era menor que la que tú encontraste, pero eso dependerá de la parametrización utilizada)
y esto para mi código 6S alterado:
- Modelo: 0.4543900000
- Fórmula: 0.4543942552
- Perc Diff: -0.0009364659%
Mucho mejor.
Como referencia, para investigar esto he utilizado Py6S la interfaz de Python para el modelo 6S que escribí. Usé las siguientes funciones para calcular automáticamente los resultados usando la fórmula de un Py6S SixS y calcular automáticamente la diferencia porcentual:
def calc_acr(radiance, xa, xb, xc):
y = xa * radiance - xb
acr = y/(1.0 + xc * y)
return acr
def calc_acr_from_obj(radiance, s):
return calc_acr(radiance, s.outputs.coef_xa, s.outputs.coef_xb, s.outputs.coef_xc)
def difference_between_formula_and_model(s):
formula = calc_acr_from_obj(s.outputs.measured_radiance, s)
model = s.outputs.atmos_corrected_reflectance_lambertian
diff = model - formula
perc_diff = (diff / model) * 100
print("Model: %.10f" % model)
print("Formula: %.10f" % formula)
print("Perc Diff: %.10f%%" % perc_diff)
y mis errores de ejemplo anteriores provienen de ejecutar Py6S utilizando la siguiente parametrización:
s = SixS()
s.altitudes.set_sensor_satellite_level()
s.atmos_corr = AtmosCorr.AtmosCorrLambertianFromRadiance(240)
s.wavelength = Wavelength(PredefinedWavelengths.LANDSAT_OLI_B1)
s.run()
Sólo como una pequeña adición, si usted está corrigiendo atmosféricamente los datos de Sentinel-2 con 6S entonces usted podría considerar el uso de ARCSI - una herramienta de corrección atmosférica que utiliza Py6S internamente, pero que hace gran parte del trabajo duro por ti. La mejor manera de aprender ARCSI es con su documento tutorial .
