GIS MCP Server

Instructions

Implementation Reference

• src/gis_mcp/pysal_functions.py:1463-1660 (handler)

  The handler function implementing the 'gm_lag' tool, which fits a spatial lag model using spreg.GM_Lag (GMM/2SLS). It handles data loading, spatial weights construction (Queen, Rook, DistanceBand), instrument generation, robust inference, and spatial diagnostics.

  @gis_mcp.tool()
  
  def gm_lag(
      shapefile_path: str,
      y_col: str,                               # dependent variable
      x_cols: Union[str, List[str]],            # exogenous regressors (no constant)
      target_crs: str = "EPSG:4326",
      weights_method: str = "queen",            # 'queen'|'rook'|'distance'
      distance_threshold: float = 100000,       # meters; auto→degrees for EPSG:4326
      # IV/GMM config
      w_lags: int = 1,                          # instruments: WX, WWX, ...
      lag_q: bool = True,                       # also lag external instruments q
      yend_cols: Union[None, str, List[str]] = None,  # other endogenous regressors (optional)
      q_cols: Union[None, str, List[str]] = None,     # their external instruments (optional)
      # Inference options
      robust: Union[None, str] = None,          # None | 'white' | 'hac'
      hac_bandwidth: float = None,              # only used if robust='hac'
      spat_diag: bool = True,                   # AK test
      sig2n_k: bool = False,                    # variance uses n-k if True
      drop_na: bool = True                      # drop rows with NA in y/x/yend/q
  ) -> Dict[str, Any]:
      """
      Run spreg.GM_Lag (spatial 2SLS / GMM-IV spatial lag model) on a cross-section.
  
      Returns coefficients with SE/z/p, fit metrics, (optional) AK test, and a small data preview.
      """
      try:
          # --- sanitize inputs ---
          for s in ("shapefile_path","target_crs","weights_method"):
              v = locals()[s]
              if isinstance(v, str):
                  locals()[s] = v.replace("`","")
  
          if isinstance(x_cols, str):
              x_cols_list = [c.strip() for c in x_cols.split(",") if c.strip()]
          else:
              x_cols_list = list(x_cols)
  
          if isinstance(yend_cols, str):
              yend_cols_list = [c.strip() for c in yend_cols.split(",") if c.strip()]
          elif yend_cols is None:
              yend_cols_list = None
          else:
              yend_cols_list = list(yend_cols)
  
          if isinstance(q_cols, str):
              q_cols_list = [c.strip() for c in q_cols.split(",") if c.strip()]
          elif q_cols is None:
              q_cols_list = None
          else:
              q_cols_list = list(q_cols)
  
          if not os.path.exists(shapefile_path):
              return {"status": "error", "message": f"Shapefile not found: {shapefile_path}"}
          if not x_cols_list:
              return {"status": "error", "message": "x_cols must include at least one regressor."}
  
          # --- load + project ---
          gdf = gpd.read_file(shapefile_path)
          needed = [y_col] + x_cols_list + (yend_cols_list or []) + (q_cols_list or [])
          missing = [c for c in needed if c not in gdf.columns]
          if missing:
              return {"status": "error", "message": f"Columns not found: {missing}"}
  
          gdf = gdf.to_crs(target_crs)
  
          # --- coerce numerics & subset ---
          gdf[needed] = gdf[needed].apply(pd.to_numeric, errors="coerce")
          data = gdf[needed + [gdf.geometry.name]].copy()
  
          if drop_na:
              before = data.shape[0]
              data = data.dropna(subset=needed)
              after = data.shape[0]
              if after == 0:
                  return {"status": "error", "message": "All rows dropped due to NA in y/x/yend/q."}
              if after < before:
                  gdf = gdf.loc[data.index].reset_index(drop=True)
                  data = data.reset_index(drop=True)
  
          # --- arrays for spreg ---
          y = data[[y_col]].to_numpy(dtype=float)              # (n,1)
          X = data[x_cols_list].to_numpy(dtype=float)          # (n,k), no constant
          YEND = None if not yend_cols_list else data[yend_cols_list].to_numpy(dtype=float)
          Q = None if not q_cols_list else data[q_cols_list].to_numpy(dtype=float)
  
          # --- spatial weights ---
          import libpysal
          wm = weights_method.lower()
          if wm == "queen":
              w = libpysal.weights.Queen.from_dataframe(gdf.loc[data.index], use_index=True)
          elif wm == "rook":
              w = libpysal.weights.Rook.from_dataframe(gdf.loc[data.index], use_index=True)
          elif wm == "distance":
              thr = distance_threshold if target_crs.upper() != "EPSG:4326" else distance_threshold / 111000.0
              w = libpysal.weights.DistanceBand.from_dataframe(
                  gdf.loc[data.index], threshold=thr, binary=False, use_index=True
              )
          else:
              return {"status":"error","message":f"Unknown weights_method: {weights_method}"}
          w.transform = "r"
  
          # Drop islands and re-align data if needed
          if w.islands:
              keep = [i for i in range(len(gdf.loc[data.index])) if i not in set(w.islands)]
              if not keep:
                  return {"status":"error","message":"All units are islands under current weights."}
              # reindex everything
              y = y[keep, :]
              X = X[keep, :]
              if YEND is not None: YEND = YEND[keep, :]
              if Q is not None: Q = Q[keep, :]
              sub_gdf = gdf.loc[data.index].iloc[keep]
              if wm == "queen":
                  w = libpysal.weights.Queen.from_dataframe(sub_gdf, use_index=True)
              elif wm == "rook":
                  w = libpysal.weights.Rook.from_dataframe(sub_gdf, use_index=True)
              else:
                  thr = distance_threshold if target_crs.upper() != "EPSG:4326" else distance_threshold / 111000.0
                  w = libpysal.weights.DistanceBand.from_dataframe(sub_gdf, threshold=thr, binary=False, use_index=True)
              w.transform = "r"
  
          # --- HAC kernel weights if requested ---
          gwk = None
          if robust == "hac":
              # Build Kernel weights on centroids; ensure ones on diagonal
              coords = np.column_stack([gdf.loc[data.index].geometry.centroid.x,
                                        gdf.loc[data.index].geometry.centroid.y])
              bw = hac_bandwidth or (np.ptp(coords[:,0]) + np.ptp(coords[:,1])) / 20.0
              gwk = libpysal.weights.Kernel(coords, bandwidth=bw, fixed=True, function="triangular", diagonal=True)
              gwk.transform = "r"
  
          # --- fit GM_Lag ---
          try:
              from spreg import GM_Lag
          except ModuleNotFoundError:
              return {"status": "error", "message": "The 'spreg' package is not installed. Install with: pip install spreg"}
  
          reg = GM_Lag(
              y, X,
              yend=YEND, q=Q,
              w=w,
              w_lags=int(w_lags),
              lag_q=bool(lag_q),
              robust=robust,               # None | 'white' | 'hac'
              gwk=gwk,                     # required if robust='hac'
              sig2n_k=bool(sig2n_k),
              spat_diag=bool(spat_diag),
              name_y=y_col,
              name_x=x_cols_list,
              name_yend=(yend_cols_list or None),
              name_q=(q_cols_list or None),
              name_w=f"{weights_method}"
          )
  
          # --- package outputs ---
          def arr(x): return None if x is None else np.asarray(x).tolist()
          def zpack(z_list):
              # z_stat is list of (z, p)
              return [{"z": float(zp[0]), "p": float(zp[1])} for zp in (z_list or [])]
  
          # tiny preview (avoid geometry dtype issues)
          preview = gdf.loc[data.index, [y_col, *x_cols_list, gdf.geometry.name]].head(5).copy()
          preview["geometry_wkt"] = preview[gdf.geometry.name].apply(lambda g: g.wkt if g is not None else None)
          preview = pd.DataFrame(preview.drop(columns=[gdf.geometry.name])).to_dict(orient="records")
  
          result = {
              "n_obs": int(reg.n),
              "k_vars": int(reg.k),  # includes constant internally
              "dependent": y_col,
              "exog": x_cols_list,
              "endog": yend_cols_list,
              "instruments": q_cols_list,
              "weights_method": weights_method.lower(),
              "spec": {"w_lags": int(w_lags), "lag_q": bool(lag_q), "robust": robust, "sig2n_k": bool(sig2n_k)},
              "betas": [float(b) for b in np.asarray(reg.betas).flatten()],
              "beta_names": (["const"] + x_cols_list + (yend_cols_list or []) + ["W_y"])[:len(np.asarray(reg.betas).flatten())],
              "std_err": arr(reg.std_err),
              "z_stats": zpack(getattr(reg, "z_stat", None)),
              "pseudo_r2": float(getattr(reg, "pr2", np.nan)),
              "pseudo_r2_reduced": float(getattr(reg, "pr2_e", np.nan)),
              "sig2": float(getattr(reg, "sig2", np.nan)),
              "ssr": float(getattr(reg, "utu", np.nan)),
              "ak_test": arr(getattr(reg, "ak_test", None)),  # [stat, p] if spat_diag=True
              "pred_y_head": [float(v) for v in np.asarray(reg.predy).flatten()[:5]],
              "data_preview": preview
          }
  
          msg = "GM_Lag estimation completed successfully"
          if wm == "distance" and target_crs.upper() == "EPSG:4326":
              msg += f" (threshold interpreted as {distance_threshold/111000.0:.6f} degrees)."
          if robust == "hac":
              msg += f" (HAC bandwidth ~ {bw:.3f})."
  
          return {"status": "success", "message": msg, "result": result}
  
      except Exception as e:
          return {"status": "error", "message": f"Failed to run GM_Lag: {str(e)}"}