This function creates a probability (stochastic) grid 
        for trajectory model data using 'overall' method (based
        on normalization by nsteps * nparticles) or 'run' method
        (based on normalization by run).
        
        probability_grid = compute_probability([myfile1.nc, myfile2.nc],
                                               bbox = [-75, 23, -60, 45],
                                               nx = 1000, ny = 1000,
                                               method = 'overall')  
        stoch.compute_probability(['/media/sf_Python/trajectories.nc'],bbox=[-148,60,-146,61],nx=500,ny=500,method='overall',parameter='settlement')

def compute_probability(trajectory_files, bbox=None, 
                        nx=None, ny=None, method='overall', parameter='location'):
    """
        This function creates a probability (stochastic) grid 
        for trajectory model data using 'overall' method (based
        on normalization by nsteps * nparticles) or 'run' method
        (based on normalization by run).
        
        probability_grid = compute_probability([myfile1.nc, myfile2.nc],
                                               bbox = [-75, 23, -60, 45],
                                               nx = 1000, ny = 1000,
                                               method = 'overall')  
        stoch.compute_probability(['/media/sf_Python/trajectories.nc'],bbox=[-148,60,-146,61],nx=500,ny=500,method='overall',parameter='settlement')
    """
    xarray = np.linspace(float(bbox[0]), float(bbox[2]), int(nx)+1) 
    yarray = np.linspace(float(bbox[1]), float(bbox[3]), int(ny)+1)
    if method=='overall':
        prob = np.zeros((ny, nx))
        for runfile in trajectory_files:
            run = netCDF4.Dataset(runfile)
            if parameter == 'location':
                lat = run.variables['lat'][:].flatten()
                lon = run.variables['lon'][:].flatten()
                column_i, row_i = [], []
                for clon, clat in zip(lon, lat):
                    column_i.append(bisect.bisect(xarray, clon))
                    row_i.append(bisect.bisect(yarray, clat))
                    try:
                        prob[row_i[-1], column_i[-1]] += 1
                    except StandardError:
                        pass
            elif parameter == 'settlement':
                for i in range(run.variables['time'].shape[0]):
                    settle_index = np.where(run.variables['settled'][i,:]==1)
                    if len(settle_index[0]) > 0:
                        lat = run.variables['lat'][i, settle_index[0]].flatten()
                        lon = run.variables['lon'][i, settle_index[0]].flatten()
                        column_i, row_i = [], []    
                        for clon, clat in zip(lon, lat):
                            column_i.append(bisect.bisect(xarray, clon))
                            row_i.append(bisect.bisect(yarray, clat))
                            try:
                                prob[row_i[-1], column_i[-1]] += 1
                            except StandardError:
                                pass
            else:
                raise ValueError("Parameter for stochastic assessment not valid")
        shape = run.variables['time'].shape
        prob = prob / (shape[0] * len(trajectory_files)) # Assumes same # of particles
                                                         # for every run, may be bad
                                                         # assumtion
    elif method=='run':
        prob = []
        for i, runfile in enumerate(trajectory_files):
            prob.append(np.zeros((ny, nx)))
            run = netCDF4.Dataset(runfile)
            if parameter == 'location':
		        lat = run.variables['lat'][:].flatten()
		        lon = run.variables['lon'][:].flatten()
            elif parameter == 'settlement':
                settle_index = np.where(run.variables['settled'][-1,:]==1)[0]
                lat = run.variables['lat'][-1, settle_index].flatten()
                lon = run.variables['lon'][-1, settle_index].flatten()
            else:
                raise ValueError("Parameter for stochastic assessment not valid")
            column_i, row_i = [], []
            for clon, clat in zip(lon, lat):
                column_i.append(bisect.bisect(xarray, clon))
                row_i.append(bisect.bisect(yarray, clat))
                try:
	                if prob[i][row_i[-1], column_i[-1]] == 0:
	                	prob[i][row_i[-1], column_i[-1]] = 1
                except StandardError:
                    pass
        prob2 = np.zeros((ny, nx))
        for run in prob:
            prob2 = run + prob2
        prob = prob2 / len(prob)
    return prob