In [ ]:

    

%pylab notebook
import glob
import matplotlib.pyplot as plt
import numpy as np
from charistools.readers import ModisTileCube
from charistools.readers import read_tile
pylab.rcParams['figure.figsize'] = (10.0, 8.0)


    

In [ ]:

    

%cd /scratch/summit/brodzik/ti_model/snow_cover/MODSCAG_GF/v09/h23v05
%ls


    

In [ ]:

    

topDir = '/scratch/summit/brodzik/ti_model/snow_cover/MODSCAG_GF/v09/'
tileID = 'h23v05'
h23v05_cube = ModisTileCube("%s/%s/MODSCAG_GF_Snow.v0.9.%s_2001.h5" % (topDir, tileID, tileID), 
                            varname='fsca')
tileID = 'h23v04'
h23v04_cube = ModisTileCube("%s/%s/MODSCAG_GF_Snow.v0.9.%s_2001.h5" % (topDir, tileID, tileID), 
                            varname='fsca')


    

In [ ]:

    

doyrange = np.arange(3)+41
doyrange


    

In [ ]:

    

fig, ax = plt.subplots(2,3)
for i, doy in enumerate(doyrange):
    upper_data = h23v04_cube.read(doy=doy)
    lower_data = h23v05_cube.read(doy=doy)
    ax[0,i].imshow(upper_data, cmap=plt.cm.gray, vmin=0.0, vmax=1.0, interpolation='None')
    ax[0,i].set_title("h23v04 day %d" % doy)
    ax[0,i].axis('off')
    ax[1,i].imshow(lower_data, cmap=plt.cm.gray, vmin=0.0, vmax=1.0, interpolation='None')
    ax[1,i].set_title("h23v05  day %d" % doy)
    ax[1,i].axis('off')
#fig.savefig('/Users/brodzik/tmp/scag_gf_2001120.v6v8.png')


    

In [ ]:

    

basinDir = '/scratch/summit/brodzik/ti_model/basins/partner_basins/'
%cd $basinDir
basins = glob.glob("BL*")
basins


    

In [ ]:

    

for i, basin in enumerate(basins):
    print(basin)
    fileList = glob.glob("%s/%s/%s*basin_mask*tif" % (basinDir, basin, basin))
    print(fileList[0])
    
    mask = read_tile(fileList[0], verbose=True)
    
    if 0 == i:
        all = mask.copy()
    else:
        idx = mask == 1
        all[idx] = i + 1
        
    print(np.amin(all), np.amax(all))


    

In [ ]:

    

doy = 41
fig, ax = plt.subplots(1,2)
scag_data = h23v04_cube.read(doy=doy)
ax[0].imshow(scag_data, cmap=plt.cm.gray, vmin=0.0, vmax=1.0, interpolation='None')
ax[0].set_title("h23v04 day %d" % doy)
ax[0].axis('off')
ax[1].imshow(all, cmap=plt.cm.terrain, vmin=0.0, vmax=5.0, interpolation='None')
ax[1].set_title("BL basins")
ax[1].axis('off')
plt.tight_layout()
#fig.savefig('/projects/brodzik/charis_ti_melt/BL_basins_vs_scag_area.png')


    

In [ ]:

    

15 * 14


    

In [ ]:

    

new_cube.close()
orig_cube.close()


    

In [ ]:

    

d5.dtype, orig_d5.dtype


    

In [ ]:

    

diff = d5 - orig_d5
print(np.amin(diff), np.amax(diff))


    

In [ ]:

    

fig, ax = plt.subplots(1)
ax.imshow(diff, cmap=plt.cm.gray, interpolation='None')
ax.set_title("Diff (new - orig)")
plt.axis('off')


    

In [ ]:

    

print(d5[2397:2399, 2397:2399])


    

In [ ]:

    

print(orig_d5[2397:2399, 2397:2399])


    

In [ ]: