4

I have multiple files in one folder which are called, for example: 049.h11v11.hdf, 049.h11v12.hadf, 049.h11v13.hdf, 065.h11v11.hdf, 065.h12v12.hdf, and so on. The first 3 numbers are the julian day of MODIS LAI (MOD15A2) product.

The objective is to do daily mosaics with about 3 years of modis MOD13A2 products. At this point I have a folder named "LAI", that contains the years folder, e.g 2000, and so on. Those folder contains the .hdf files I named in the last paragraph. So, the idea is that the python code can run for all folders and get all same days and do the mosaic.

So far I've though about something like this:

#Libraries
import gdal
from gdalconst import *
import osgeo
import osr
import osgeo.gdal as gdal
import osgeo.osr as osr
import numpy as np
import glob
import os

#Roots:
RutaLai="/mnt/discoenred/jorge/scripts/MOD15_A2_mosaic/LAI/"

#Create dir
if not os.path.exists(RutaLai): os.mkdir(RutaLai)
.
.
#in between this are the reprojection and the multiplication by the LAI factor for each .hdf file.
.
.
#Here is were the mosaic code starts:
for raiz, carpetas, archivos in os.walk(RutaLai, topdown=False):

    arreglo = carpetas

diaTemporal = None  
diaVariable = []
diaArreglo = []

for carpeta in arreglo: 
    os.chdir(RutaLai+carpeta) #carpeta is defined as as the year taken from the MOD13A2 name, i.e carpeta = x[9:13]+"/" ; were x goes inside a for x in lista2, that just take the year inside the name of the file. (it is inside the code for the LAI factor, I don't put it because its really huge. If it is really needed then I can put it)
    #Lai images lists
    lista3=glob.glob("*.hdf")
    lista3.sort()

    for archivo in lista3:
        dia = archivo.split(".")

        #Definition of temporal day (diaTemporal) variable to avoid error.
        if (diaTemporal == None):
            diaTemporal = dia[0]

        if (diaTemporal == dia[0]):
            diaVariable.append( dia[0] )

        else:
            if (diaVariable != None):
                diaArreglo.append( diaVariable )

        diaTemporal = dia[0]

    #RutaDestinoMosaico = RutaMosaic + carpeta + "mosaico_" + ".hdf"    
    #os.system("gdal_merge.py -o "+ RutaDestinoMosaico + " -of ENVI -n 0 '"+archivo1+"' '"+archivo2+"'")

The idea is that the code go to all e.g 049 files in one folder and do the mosaic that will be saved in a folder inside the year called "mosaic".. the files should go like "mosaic_049.hdf" (for the example given above). In order to understand it better I took a screen shot of the directory:

LAI files with the scale factor

It is important to notice that the picture (inside year "2000") contains the files that are inputs for the mosaic. For this all I'm using a Linux machine with Ubuntu.

2

So you just want to mosaic all the tiles for a given day? That's a perfect Job for GDALs VRTs.

gdalbuildvrt mosaic_049.vrt 049*.hdf

or from Python for all days

import subprocess
import glob

for day in range(0, 365):
    day = str('%0.3d' % day)
    cmd = ["gdalbuildvrt", "mosaic_"+day+".vrt", glob.glob(day+"*.hdf"]
    subprocess.call(cmd)

edit: included cleaner day formatting by @Detlev

edit2: If you are using GDAL >2.0 you can call buildvrt as a function.

from osgeo import gdal
import glob


for day in range(0, 365):
    day = str('%0.3d' % day)
    gdal.BuildVRT("mosaic_"+day+".vrt", glob.glob(day+".hdf")
  • 1
    Abbreviate lines 5-9 to: day = str('%0.3d' % day) – Detlev Jan 28 '16 at 15:23
  • @Detlev great hint. I knew I was doing something very unpythonic when I used two consecutive if statements - just couldn't figure out the string formatting fast enough. – Kersten Jan 28 '16 at 15:30
  • @Kersten, sorry for not answering before. Thank you for your time, but this did not solve the problem. First of all, I need to make it with gdal merge. I know I can use then gdal translate to convert de .vrt file into .hdr, but the thing is that I'm doing this for more than 20.000 images, so I really need to optimize code. I tried changing gdal buildvrt for gdal merge with no exit. Anyway, thank you, the day in range seems a pretty good idea, hope it works. – Jorge V Feb 3 '16 at 16:56
  • 1
    When you want to merge the files inside the VRT into a .hdr it would be as simple as calling gdal_translate -of ENVI mosaic_049.vrt mosaic_049.hdr. gdal_merge uses VRTs internally anyway if I remember correctly so I don't think you'll see a performance increase in calling it directly. – Kersten Feb 3 '16 at 19:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.