from scipy.interpolate import interp1d
rc('font', **{'family': 'serif', 'serif': ['Computer Modern']})
rc('font', size=20)
rc('text', usetex=True)
rc('xtick', labelsize=20, color='black')
rc('ytick', labelsize=20, color='black')
rc('lines', linewidth=3.0)
rcParams['axes.linewidth'] = 1.0
rcParams['axes.edgecolor'] = 'black'
rcParams['figure.figsize'] = (10.0, 8.0)
rc('text', usetex=True)


xyd = loadtxt("hypsometric_curve.out")
xy  = (zeros(99),zeros(99)) 
xy[0][:] = xyd[:,0]; xy[1][:] = xyd[:,1]

def hypsometric_curve(level):
  data  = loadtxt("etopo1_b_c");
  len_longi = 360*6+1;
  len_lati  = 180*6+1;
  longi = data[:,0].reshape(len_lati,len_longi);
  lati  = data[:,1].reshape(len_lati,len_longi);
  elev  = data[:,2].reshape(len_lati,len_longi);
  area_grid     = zeros((len_lati,len_longi))
  elev_grid     = zeros((len_lati,len_longi))
  dlat  = (lati[1,0]-lati[0,0])  *pi/180.   # in radians
  dlong = (longi[0,1]-longi[0,0])*pi/180.   # in radians
  R     = 6378137. # Radius of the earth
  for i in xrange(len_lati):
    for j in xrange(len_longi):
      lati_ij   = (lati[i,j] ) *pi/180. # in radians  
      longi_ij  = (longi[i,j]) *pi/180. # in radians  
      dy        = R*dlat
      dx        = R*cos(lati_ij)*dlong
      area_grid[i,j]   = fabs(dx*dy)  # Area of the grid
      elev_grid[i,j]   = elev[i,j]

  area_grid = area_grid.reshape((len_longi)*(len_lati))
  area_grid = area_grid/sum(area_grid)
  elev_grid = elev_grid.reshape((len_longi)*(len_lati))

  # Calculate the histogram
  binsl = linspace(min(elev_grid), max(elev_grid), level)
  hist , bin_edges= histogram(elev_grid, bins=binsl, weights=area_grid, normed=True, density=True)

  # calculate the hypsometric curve
  xp = zeros(level-1)
  yp = zeros(level-1)
  for i in xrange(level-1):
    xp[i] = sum(hist[level-2-i:level-1]*diff(bin_edges[level-i-2:level]))
    yp[i] = (bin_edges[level-1-i]+bin_edges[level-2-i])/2.
    
  return (xp, yp)

def surface_to_elev (point, x, y):
    f_value = interp1d(x, y)
    return f_value(point)
    

def elev_to_surface (point, x, y):
    from scipy.interpolate import interp1d
    f_value = interp1d(y, x)
    return f_value(point)
    

def plot_hypsometric_curve(x,y):
    plot(x,y,'k')
    x0 = elev_to_surface(0.,x,y)
    still_water = linspace(x0,max(x),100)
    bathy       = surface_to_elev(still_water,x,y)
    fill_between(x, min(y), y, color='darkgoldenrod')
    fill_between(still_water, 0, bathy, color='b')
    xlim(min(x),max(x))
    ylim(min(y),max(y))
    xlabel(r'Surface/Total Surface')
    ylabel(r'Elevation')

def separate_land_ocean(elv,per):
    nx = len(elv)
    new_ele_oc=[]
    new_perc_oc=[]
    new_ele_la=[]
    new_perc_la=[]
    for i in range(nx):
        if elv[i]<=-0.0:
            new_ele_oc.append(elv[i])
            new_perc_oc.append(per[i])
        if elv[i]>-0.0:
            new_ele_la.append(elv[i])
            new_perc_la.append(per[i])
    new_ele_oc=array(new_ele_oc)
    new_perc_oc=array(new_perc_oc)
    new_ele_la=array(new_ele_la)
    new_perc_la=array(new_perc_la)
    return new_ele_oc, new_perc_oc, new_ele_la, new_perc_la