#converts hdf files to tiffs, 

#slicing is the number of pixels that will be kept from original input

#slicing extent is: [xmin,xmax,ymin,ymax], if not altered this will automatically cover the entire raster



#CAREFUL: lines 38-41 are input data specific and may need to be changed manually

#line 38 changes pixel values (eg if a conversion is necessary), the other three are raster pixel size and extent



import os, gdal

from gdalconst import *



def hdfTOtif(nameHDF,outFile, slicing = [0,0,0,0], projName = "EPSG:4326"):  # Standard coordinate system is assumed to be WGS84



    driver = gdal.GetDriverByName('hdf4')

    driver2 = gdal.GetDriverByName('Gtiff')

    driver.Register()

    driver2.Register()

    

    inHDF = gdal.Open(nameHDF, GA_ReadOnly)

    cols = inHDF.RasterXSize

    rows = inHDF.RasterYSize

    

    array = inHDF.ReadAsArray(0, 0, cols, rows) #HDF to numpyArray

    

    #Define pixel ranges to slice original array 

    if slicing[0] == 0 and slicing[1] == 0 and slicing[2] == 0 and slicing[3] == 0:

        x_Start = 0

        y_Start = 0

        x_End = rows

        y_End = cols

    else:

        x_Start = slicing[0]

        x_End = slicing[1]        

        y_Start = slicing[2]

        y_End = slicing[3]

    

    array = array[y_Start:y_End,x_Start:x_End]

    

    array = array * 0.004 - 0.08 #DN values to NDVI     

    PIXEL_SIZE = 0.0089285714

    x_min = 91.9910714286 + x_Start * PIXEL_SIZE #if sliced, needs be updated

    y_max = 29.0089285714 - y_Start * PIXEL_SIZE    

    

    proj = osr.SpatialReference()

    proj.SetWellKnownGeogCS( projName ) #Get the long coordinate system name

    wkt_projection = proj.ExportToWkt()

    

    x_pixels = array.shape[1]

    y_pixels = array.shape[0]

   

    outDataset = driver2.Create(

            outFile,

            x_pixels,

            y_pixels,

            1,

            gdal.GDT_Float32, )

    

    outDataset.SetGeoTransform((

            x_min,    # 0

            PIXEL_SIZE,  # 1

            0,                      # 2

            y_max,    # 3

            0,                      # 4

            -PIXEL_SIZE))  

    

    outDataset.SetProjection(wkt_projection)

    outDataset.GetRasterBand(1).WriteArray(array)

    outDataset.FlushCache()  # Write to disk.



 

    

# This is how you use the hdfTOtif function to convert an entire folder of hdf files



inFol = ".../myInputFolder/"

outFol = ".../myOutputFolder/"



for filex in os.listdir(inFol):

    extent = [800,2300,300,2500]

    nameHDF = inFol + filex

    outFile = outFol + filex[:-3] + "tif"

    hdfTOtif(nameHDF,outFile,extent)