docs: def load_training_vector(response_shapes, explanatory_rasters, response_field, metric='mean'):
      """
      Parameters
      ----------
      response_shapes : Source of vector features for raster_stats;
                        can be OGR file path or iterable of geojson-like features
      response_field : Field name containing the known response category (must be numeric)
      explanatory_rasters : List of Paths to GDAL rasters containing explanatory variables
      metric : Statistic to aggregate explanatory data across line and polygon vector features
               Defaults to 'mean' (optional)
  
      Returns
      -------
      train_xs : Array of explanatory variables
      train_y : 1xN array of known responses
      """
      from rasterstats import zonal_stats
      all_means = []
      all_zones = None
  
      for i, raster in enumerate(explanatory_rasters):
          logger.debug("Rasters stats on %s" % raster)
  
          stats = zonal_stats(response_shapes, raster, stats=metric, prefix="pyimpute_", geojson_out=True)
  
          zones = [x['properties'][response_field] for x in stats]
          if all_zones:
              assert zones == all_zones
          else:
              all_zones = zones
  
          means = [x['properties']['pyimpute_' + metric] for x in stats]
          all_means.append(means)
  
      train_y = np.array(all_zones)
      train_xs = np.array(all_means).T
  
      return train_xs, train_y