docs: def load_training_rasters(response_raster, explanatory_rasters, selected=None):
      """
      Parameters
      ----------
      response_raster : Path to GDAL raster containing responses
      explanatory_rasters : List of Paths to GDAL rasters containing explanatory variables
  
      Returns
      -------
      train_xs : Array of explanatory variables
      train_ys : 1xN array of known responses
      """
  
      with rasterio.open(response_raster) as src:
          response_data = src.read().flatten()
  
      if selected is None:
          train_y = response_data
      else:
          train_y = response_data[selected]
  
      selected_data = []
      for rast in explanatory_rasters:
          with rasterio.open(rast) as src:
              explanatory_data = src.read().flatten()
          assert explanatory_data.size == response_data.size
          if selected is None:
              selected_data.append(explanatory_data)
          else:
              selected_data.append(explanatory_data[selected])
  
      train_xs = np.asarray(selected_data).T
      return train_xs, train_y