GIS MCP Server

Instructions

Implementation Reference

• src/gis_mcp/pysal_functions.py:882-1035 (handler)

  Handler function implementing the OLS regression tool with spatial diagnostics using libpysal and spreg libraries. Includes data loading, weight construction, model fitting, and spatial diagnostics (Moran's I on residuals). Registered via @gis_mcp.tool() decorator.

  @gis_mcp.tool()
  def ols_with_spatial_diagnostics_safe(
      data_path: str,
      y_field: str,
      x_fields: List[str],
      weights_path: Optional[str] = None,
      weights_method: str = "queen",
      id_field: Optional[str] = None,
      threshold: Optional[float] = None,
      k: Optional[int] = None,
      binary: bool = True
  ) -> Dict[str, Any]:
      """
      Safe MCP pipeline: Read shapefile, build/load W, convert numeric, check NaNs, run OLS.
  
      Parameters:
      - data_path: path to shapefile or GeoPackage
      - y_field: dependent variable column name
      - x_fields: list of independent variable column names
      - weights_path: optional path to existing weights file (.gal or .gwt)
      - weights_method: 'queen', 'rook', 'distance_band', or 'knn' (used if weights_path not provided)
      - id_field: optional attribute name to use as observation IDs
      - threshold: required if method='distance_band'
      - k: required if method='knn'
      - binary: True for binary weights (DistanceBand only)
      """
      try:
          # --- Step 1: Read data ---
          if not os.path.exists(data_path):
              return {"status": "error", "message": f"Data file not found: {data_path}"}
  
          gdf = gpd.read_file(data_path)
          if gdf.empty:
              return {"status": "error", "message": "Input file contains no features"}
  
          # --- Step 2: Extract and convert y and X ---
          if y_field not in gdf.columns:
              return {"status": "error", "message": f"Dependent variable '{y_field}' not found in dataset"}
          if any(xf not in gdf.columns for xf in x_fields):
              return {"status": "error", "message": f"Independent variable(s) {x_fields} not found in dataset"}
  
          y = gdf[[y_field]].astype(float).values
          X = gdf[x_fields].astype(float).values
  
          # --- Step 3: Check for NaNs or infinite values ---
          if not np.all(np.isfinite(y)):
              return {"status": "error", "message": "Dependent variable contains NaN or infinite values"}
          if not np.all(np.isfinite(X)):
              return {"status": "error", "message": "Independent variables contain NaN or infinite values"}
  
          # --- Step 4: Load or build weights ---
          import libpysal
          if weights_path:
              if not os.path.exists(weights_path):
                  return {"status": "error", "message": f"Weights file not found: {weights_path}"}
              w = libpysal.open(weights_path).read()
          else:
              coords = [(geom.x, geom.y) for geom in gdf.geometry]
              wm = weights_method.lower()
              if wm == "queen":
                  w = libpysal.weights.Queen.from_dataframe(gdf, idVariable=id_field)
              elif wm == "rook":
                  w = libpysal.weights.Rook.from_dataframe(gdf, idVariable=id_field)
              elif wm == "distance_band":
                  if threshold is None:
                      return {"status": "error", "message": "Threshold is required for distance_band"}
                  ids = gdf[id_field].tolist() if id_field and id_field in gdf.columns else None
                  w = libpysal.weights.DistanceBand(coords, threshold=threshold, binary=binary, ids=ids)
              elif wm == "knn":
                  if k is None:
                      return {"status": "error", "message": "k is required for knn"}
                  ids = gdf[id_field].tolist() if id_field and id_field in gdf.columns else None
                  w = libpysal.weights.KNN(coords, k=k, ids=ids)
              else:
                  return {"status": "error", "message": f"Unsupported weights method: {weights_method}"}
  
          w.transform = "r"  # Row-standardize for regression
  
          # --- Step 5: Fit OLS with spatial diagnostics ---
          try:
              from spreg import OLS
          except ModuleNotFoundError:
              return {"status": "error", "message": "The 'spreg' package is not installed. Install with: pip install spreg"}
          
          # Ensure y is shape (n, 1) and X is (n, k) for spreg
          y = y.reshape(-1, 1) if len(y.shape) == 1 else y
          if X.shape[1] != len(x_fields):
              X = X.T if X.shape[0] == len(x_fields) else X
          
          # Run OLS with spatial diagnostics
          ols_model = OLS(y, X, w=w, spat_diag=True, name_y=y_field, name_x=x_fields, name_ds="dataset")
  
          # --- Step 6: Collect results ---
          # Extract Moran's I for residuals if available (from diagnostics)
          moran_residual = None
          moran_pvalue = None
          if hasattr(ols_model, "diagnostics") and ols_model.diagnostics is not None:
              diag = ols_model.diagnostics
              # Check for Moran's I in diagnostics
              if isinstance(diag, dict):
                  if "moran_res" in diag and diag["moran_res"] is not None:
                      mor_res = diag["moran_res"]
                      if len(mor_res) >= 2:
                          moran_residual = float(mor_res[0])
                          moran_pvalue = float(mor_res[1])
          
          # Extract beta names
          beta_names = []
          if hasattr(ols_model, "name_x") and ols_model.name_x:
              beta_names = list(ols_model.name_x)
          else:
              beta_names = [f"x_{i}" for i in range(len(x_fields))]
          
          # Add constant to names if present
          has_constant = hasattr(ols_model, "constant") and ols_model.constant
          if has_constant:
              beta_names = ["constant"] + beta_names
          
          # Extract betas
          betas_dict = {}
          if hasattr(ols_model, "betas") and ols_model.betas is not None:
              betas_flat = ols_model.betas.flatten()
              # Match betas to names
              for i, beta_val in enumerate(betas_flat):
                  if i < len(beta_names):
                      betas_dict[beta_names[i]] = float(beta_val)
                  else:
                      betas_dict[f"beta_{i}"] = float(beta_val)
          
          # Extract standard errors
          std_err_list = []
          if hasattr(ols_model, "std_err") and ols_model.std_err is not None:
              std_err_list = ols_model.std_err.tolist() if hasattr(ols_model.std_err, "tolist") else list(ols_model.std_err)
          
          results = {
              "n_obs": int(ols_model.n) if hasattr(ols_model, "n") else int(len(y)),
              "r2": float(ols_model.r2) if hasattr(ols_model, "r2") else None,
              "std_error": std_err_list,
              "betas": betas_dict,
              "moran_residual": moran_residual,
              "moran_pvalue": moran_pvalue,
          }
  
          return {
              "status": "success",
              "message": "OLS regression with spatial diagnostics completed successfully",
              "result": results,
              "regression_results": results  # Also include as regression_results for test compatibility
          }
  
      except Exception as e:
          logger.error(f"Error in ols_with_spatial_diagnostics_safe: {str(e)}")
          return {"status": "error", "message": f"Failed to run OLS regression: {str(e)}"}

• src/gis_mcp/pysal_functions.py:882-882 (registration)

  The @gis_mcp.tool() decorator registers this function as an MCP tool named 'ols_with_spatial_diagnostics_safe'. The FastMCP instance 'gis_mcp' is defined in src/gis_mcp/mcp.py.

  @gis_mcp.tool()

• src/gis_mcp/pysal_functions.py:108-134 (helper)

  Helper function for loading data and creating spatial weights, used in several esda tools but not directly in this OLS tool.

  def pysal_load_data(shapefile_path: str, dependent_var: str, target_crs: str, distance_threshold: float):
      """Common loader and weight creation for esda statistics."""
      if not os.path.exists(shapefile_path):
          return None, None, None, None, f"Shapefile not found: {shapefile_path}"
  
      gdf = gpd.read_file(shapefile_path)
      if dependent_var not in gdf.columns:
          return None, None, None, None, f"Dependent variable '{dependent_var}' not found in shapefile columns"
  
      gdf = gdf.to_crs(target_crs)
  
      effective_threshold = distance_threshold
      unit = "meters"
      if target_crs.upper() == "EPSG:4326":
          effective_threshold = distance_threshold / 111000
          unit = "degrees"
  
      y = gdf[dependent_var].values.astype(np.float64)
      import libpysal
      w = libpysal.weights.DistanceBand.from_dataframe(gdf, threshold=effective_threshold, binary=False)
      w.transform = 'r'
  
      for island in w.islands:
          w.weights[island] = [0] * len(w.weights[island])
          w.cardinalities[island] = 0
  
      return gdf, y, w, (effective_threshold, unit), None