Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
A script to calculate areas of perennial snow from MOD10A2 data
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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.