A script to calculate areas of perennial snow from MOD10A2 data

perrenialSnow.py

import os, subprocess, glob, numpy
from datetime import datetime
import gdal_calculations as gdalcalc

def hdf_dir_to_tif(input_dir):
    acquisition_date = input_dir[-11:-1]
    print 'Converting HDFs to TIFs'
    for hdf_filepath in glob.glob(input_dir+'*.hdf'):
        tif_filepath = hdf_filepath +'.tif'
        subprocess.call(
            'gdal_translate -co COMPRESS=PACKBITS HDF4_EOS:EOS_GRID:{0}:MOD_Grid_Snow_500m:Maximum_Snow_Extent {1}'.format(hdf_filepath, tif_filepath),
            shell = True
            )
        os.remove(hdf_filepath)
    for hdf_filepath in glob.glob(input_dir+'*.jpg'):
        os.remove(hdf_filepath)
    for hdf_filepath in glob.glob(input_dir+'*.xml'):
        os.remove(hdf_filepath)
    print 'Merging TIFs'
    subprocess.call('gdal_merge.py -co COMPRESS=PACKBITS -o {0}{1}.tif {0}*hdf.tif'.format(input_dir, acquisition_date), shell=True)
    for hdf_filepath in glob.glob(input_dir+'*hdf.tif'):
        os.remove(hdf_filepath)
        
def day_files_to_sum_files(input_datasets, testvals):
    """
    input_datasets should be a list of tif file paths.
    This function adds up the 200 values by pixel across the datasets, and also adds up the 25 values.
    The output Dataset has the 200 values in band 1 and 25 in band 2.
    """
    # load tifs
    tif_data={}
    for tif_filepath in input_datasets:
        acquisition_date = datetime.strptime(tif_filepath.split('/')[-1],('%Y.%m.%d.tif'))
        tif_data[acquisition_date]=gdalcalc.Dataset(tif_filepath)
        print "loaded data for {0}".format(acquisition_date)

    #test values
    for testval in testvals:
        summer = []
        winter = []
        for acquisition_date in tif_data:
            print "finding pixels for " + str(acquisition_date) + " " + str(testval)
            tested = tif_data[acquisition_date] == testval
            if 3 < acquisition_date.month < 9:
                summer.append(tested)
            else:
                winter.append(tested)

        print 'there are {0} summer datasets'.format(len(summer))
        print 'calculating summer sum for {0}'.format(testval)
        loop_inplace_sum(summer).save('summer_sum_of_{0}.tif'.format(testval))

        print 'there are {0} winter datasets'.format(len(winter))
        print 'calculating winter sum for {0}'.format(testval)
        loop_inplace_sum(summer).save('winter_sum_of_{0}.tif'.format(testval))
        
def persnow_from_dir(input_dir):
    layers = []
    for val in [200, 25]:
        for season in ['summer', 'winter']:
           layers.append(gdalcalc.Dataset('{0}_sum_of_{1}.tif'.format(season, val)))
    return [persnow(*layers), undetermined_for_persnow(*layers)]

def persnow(summer_snow_sum, winter_snow_sum, summer_nosnow_sum, winter_nosnow_sum):
    persnow = summer_snow_sum & winter_snow_sum & (summer_nosnow_sum == 0) & (winter_nosnow_sum == 0)
    return persnow.save('persnow.tif')

def undetermined_for_persnow(summer_snow_sum, winter_snow_sum, summer_nosnow_sum, winter_nosnow_sum):
    undet = (summer_nosnow_sum == 0) & (winter_nosnow_sum == 0) & (winter_snow_sum == 0) & (winter_nosnow_sum == 0)
    return undet.save('persnow_undet.tif')

def loop_inplace_sum(arrlist):
    # assumes len(arrlist) > 0
    sum = arrlist[0].copy()
    for a in arrlist[1:]:
        sum += a
    return sum