In [ ]:
%pylab notebook
import csv
import datetime as dt
import glob
import matplotlib.pyplot as plt
import matplotlib.dates as md
from nose.tools import set_trace
from charistools.convertors import Dem2Hypsometry
from charistools.hypsometry import Hypsometry
from charistools.meltModels import TriSurfTempIndexMelt
from charistools.meltModels import ImshowTriSurfMelt
from charistools.meltModels import PlotTriSurfInput
from charistools.meltModels import PlotTriSurfMelt
from charistools.modelEnv import ModelEnv
from charistools.timeSeries import TimeSeries
import pandas as pd
import re
import os
In [ ]:
majorConfigFile = '/Users/brodzik/ipython_notebooks/charis/major_basins_modelEnv_config.ini'
majorBasinEnv = ModelEnv(tileConfigFile=majorConfigFile)
snowyConfigFile = '/Users/brodzik/ipython_notebooks/charis/snowy_basins_modelEnv_config.ini'
snowyBasinEnv = ModelEnv(tileConfigFile=snowyConfigFile)
In [ ]:
snowyBasinEnv.tileConfig
In [ ]:
# ids = ['AM', 'BR', 'GA_v01', 'IN_v01', 'SY_v01']
ids = ['SY_v01']
snowyBasinIDs = ["%s_snow_fullbasinmasks" % id for id in ids]
fullBasinIDs = ["%s_fullbasinmasks" % id for id in ids]
In [ ]:
ids, fullBasinIDs, snowyBasinIDs
In [ ]:
#help(Hypsometry)
In [ ]:
# Generate a list of all OBJECTID basins
def objectIDs(id):
%cd /Users/brodzik/projects/CHARIS/derived_hypsometries/MODSCAG_GF_v09_fromFile_rainfall_less_ET/
objectIDs = []
ids = glob.glob("%s_OBJECTID*" % id)
ids.sort()
print("There are %d OBJECTID sub-basins for the %s major basin" % (len(ids), id))
return(ids)
In [ ]:
ids
In [ ]:
years = np.arange(2) + 2008
years
In [ ]:
for year in years:
for i, id in enumerate(ids):
# Read the full basin elevation and total the area
fullFile = majorBasinEnv.hypsometry_filename(type='area_by_elevation', drainageID=fullBasinIDs[i])
full_hyps = Hypsometry(fullFile)
full_basin_area = full_hyps.data_by_doy().loc['NoDate']
print("%s_area_km2: %.4f" % (fullBasinIDs[i], full_basin_area))
# Read the full basin elevation and total the area
snowyFile = majorBasinEnv.hypsometry_filename(type='area_by_elevation', drainageID=snowyBasinIDs[i])
snowy_hyps = Hypsometry(snowyFile)
snowy_basin_area = snowy_hyps.data_by_doy().loc['NoDate']
print("%s_area_km2: %.4f" % (snowyBasinIDs[i], snowy_basin_area))
# Get a list of the snowy_basins for this id
# Read each area file for each objectID in turn, and sub-total by doy
objects = objectIDs(id)
#objects = objects[:3]
partition_types = ['snow_on_land', 'exposed_glacier_ice', 'snow_on_ice']
sol = pd.DataFrame()
soi = pd.DataFrame()
egi = pd.DataFrame()
for object in objects:
print("Next objectID: %s" % object)
SOLFile = snowyBasinEnv.hypsometry_filename(
type='snow_on_land_by_elevation',
drainageID=object,
year=year,
modice_nstrikes=3)
SOIFile = snowyBasinEnv.hypsometry_filename(
type='snow_on_ice_by_elevation',
drainageID=object,
year=year,
modice_nstrikes=3,
threshold='fromFile',
ablation_method='grsize_scag')
EGIFile = snowyBasinEnv.hypsometry_filename(
type='exposed_glacier_ice_by_elevation',
drainageID=object,
year=year,
modice_nstrikes=3,
threshold='fromFile',
ablation_method='grsize_scag')
solHyps = Hypsometry(SOLFile)
soiHyps = Hypsometry(SOIFile)
egiHyps = Hypsometry(EGIFile)
sol_by_doy = solHyps.data_by_doy()
sol_by_doy.name = object
soi_by_doy = soiHyps.data_by_doy()
soi_by_doy.name = object
egi_by_doy = egiHyps.data_by_doy()
egi_by_doy.name = object
sol = sol.append(sol_by_doy)
soi = soi.append(soi_by_doy)
egi = egi.append(egi_by_doy)
# Make a new DataFrame with the areas by date and SOL/SOI/EGI columns
df = pd.concat([sol.sum(), egi.sum(), soi.sum()], axis=1)
df.columns = ['SOL_km2', 'EGI_km2', 'SOI_km2']
df['Snow_ice_area_km2'] = df.sum(axis=1)
df['Full_basin_area_km2'] = full_basin_area
df['Full_basin_less_snow_ice_area_km2'] = full_basin_area - df['Snow_ice_area_km2']
df['Snowy_basin_area_km2'] = snowy_basin_area
df['Snowy_basin_less_snow_ice_area_km2'] = snowy_basin_area - df['Snow_ice_area_km2']
# Write out the data frame
outFile = os.path.join(
"/Users/brodzik/projects/CHARIS/derived_hypsometries/MODSCAG_GF_v09_fromFile_rainfall_less_ET",
"%s.%4d.partition_area.new.txt" % (id, year))
print("Outfile will be: %s" % outFile)
fh = open(outFile, 'w')
fh.write("# Major_basin: %s\n" % id)
fh.write("# Snowy basins: %s\n" % (", ".join(objects)))
fh.write("# Number_of_snowy_basins: %d\n" % (len(objects)))
fh.write(str(len(df.columns)) + '\n')
fh.write(' '.join([str(col) for col in df.columns]) + '\n')
fh.close()
df['yyyy'] = df.index.year
df['mm'] = df.index.month
df['dd'] = df.index.day
df['doy'] = df.index.dayofyear
df = df.reindex_axis(
['yyyy', 'mm', 'dd', 'doy'] +
list(df.columns[:-4]), axis=1)
format = "%.4f"
df.to_csv(outFile, mode='a', header=False, index=False, sep=" ",
float_format=format, quoting=csv.QUOTE_NONE)
In [ ]:
%cat /Users/brodzik/projects/CHARIS/derived_hypsometries/MODSCAG_GF_v09_fromFile_rainfall_less_ET/SY_v01.2001.partition_area.txt
In [ ]:
fig, axes = plt.subplots(1)
for id in majorBasinIDs:
outFile = myEnv.hypsometry_filename(type='area_by_elevation', drainageID=id)
print(outFile)
hyps = Hypsometry(outFile)
hyps.data.loc['NoDate'].plot(label=id)
axes.legend()
In [ ]:
fig, axes = plt.subplots(1)
for id in majorBasinIDs:
outFile = myEnv.hypsometry_filename(type='area_by_elevation', drainageID=id)
print(outFile)
hyps = Hypsometry(outFile)
hyps.data.loc['NoDate'].plot(label=id)
axes.legend()
In [ ]:
#hyps.data
In [ ]: