docs: 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