Transformar día juliano en día calendario mediante Python

La convención para el nombre de los archivos GeoTIFF de las misiones Landsat (4, 5, 7) sigue el formato que se encuentra en este link; donde el llamado día juliano forma parte de sus características. Por tanto, en el siguiente nombre de ejemplo:

    LT50380322011235PAC01_B3.TIF

la parte resaltada en negrilla hace referencia al día 235 (día juliano) del año 2011 para este producto TM Landsat 5.

La corrección atmósférica de las imágenes, como por ejemplo la transformación a reflectancia aparente de los valores de radiacia, necesita como parámetro la distancia Tierra-Sol (expresada en unidades astronómicas) la cual, en su cálculo, exige el día específico del año en la que fue tomada la imagen. El paquete Landsat del lenguaje R requiere, para su operación, que el formato de la fecha sea el siguiente: yyyy-mm-dd.

Teniendo presente lo anterior, se desarrolló un código que transforma el conjunto año-día juliano en una fecha según el formato requerido; tomando también en consideración los años bisiestos. Su ejecución en el código siguiente se hizo para el día del ejemplo (2011235), es decir, 23 de agosto de 2011. No obstante, se probó en un número exhaustivo de situaciones y los resultados siempre fueron satistfactorios. A los interesados se les da la oportunidad de generalizarlo.

def leap_year(year):
    return  ( year % 4 == 0 and year % 100 !=0 ) or year % 400 == 0 

year = 2011
leap = leap_year(year)

days = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]

if leap == True:
    days[1] = 29

day = 235
sum = 0

for k, item in enumerate(days):
    sum += item
    
    if sum >= day:
        dd = item - (sum - day)
        
        if k+1 < 10 or dd < 10:
            mm = str(k+1).zfill(2)
            dd = str(dd).zfill(2)
        else:
            mm = str(k+1)
            dd = str(dd)

        formatted_day = str(year) + "-" + mm + "-" +  dd
        
        break

print formatted_day

Ejecutado en la Python Console de QGIS se tiene lo siguiente:

Python Console 
Use iface to access QGIS API interface or Type help(iface) for more info
execfile(u'C:/pyqgis_scripts/dia_calendario.py'.encode('mbcs'))
2011-08-23

cuya salida en R produce el resultado siguiente:

> landsat::ESdist("2011-08-23")
[1] 1.011209
Esta entrada fue publicada en Código Python. Guarda el enlace permanente.

Responder

Por favor, inicia sesión con uno de estos métodos para publicar tu comentario:

Logo de WordPress.com

Estás comentando usando tu cuenta de WordPress.com. Cerrar sesión / Cambiar )

Imagen de Twitter

Estás comentando usando tu cuenta de Twitter. Cerrar sesión / Cambiar )

Foto de Facebook

Estás comentando usando tu cuenta de Facebook. Cerrar sesión / Cambiar )

Google+ photo

Estás comentando usando tu cuenta de Google+. Cerrar sesión / Cambiar )

Conectando a %s