GIS MCP Server

Instructions

Implementation Reference

• src/gis_mcp/pysal_functions.py:1143-1321 (handler)

  The core handler function implementing the 'spatial_markov' MCP tool. It loads panel data from a shapefile, constructs spatial weights, discretizes values into Markov classes, computes conditional transition probabilities, steady states, and performs inference tests (chi2, Q, LR).
  def spatial_markov( shapefile_path: str, value_columns: Union[str, List[str]], # time-ordered, oldest -> newest target_crs: str = "EPSG:4326", weights_method: str = "queen", # 'queen'|'rook'|'distance' distance_threshold: float = 100000, # meters (converted to degrees if 4326) k: int = 5, # classes for y (quantile bins if continuous) m: int = 5, # classes for spatial lags fixed: bool = True, # pooled quantiles across all periods permutations: int = 0, # >0 to get randomization p-values for x2 relative: bool = True, # divide each period by its mean drop_na: bool = True, # drop features with any NA across time columns fill_empty_classes: bool = True # handle empty bins by making them self-absorbing ) -> Dict[str, Any]: """Run giddy Spatial Markov on a panel (n regions x t periods) from a shapefile.""" try: # --- sanitize --- for s in ("shapefile_path","target_crs","weights_method"): v = locals()[s] if isinstance(v, str): locals()[s] = v.replace("`","") if isinstance(value_columns, str): value_cols = [c.strip() for c in value_columns.replace("`","").split(",") if c.strip()] else: value_cols = list(value_columns) if not os.path.exists(shapefile_path): return {"status": "error", "message": f"Shapefile not found: {shapefile_path}"} if len(value_cols) < 2: return {"status": "error", "message": "value_columns must include at least 2 time steps (wide format)."} # --- load + project --- gdf = gpd.read_file(shapefile_path) missing = [c for c in value_cols if c not in gdf.columns] if missing: return {"status": "error", "message": f"Columns not found: {missing}"} gdf = gdf.to_crs(target_crs) # Ensure numeric gdf[value_cols] = gdf[value_cols].apply(pd.to_numeric, errors="coerce") # --- prepare Y (n x t) --- Y = gdf[value_cols].to_numpy(copy=True).astype(float) # shape (n, t) if drop_na: mask = ~np.any(np.isnan(Y), axis=1) if mask.sum() < Y.shape[0]: gdf = gdf.loc[mask].reset_index(drop=True) Y = Y[mask, :] # optional relative values (period-wise) if relative: col_means = Y.mean(axis=0) col_means[col_means == 0] = 1.0 Y = Y / col_means # --- spatial weights --- import libpysal wm = weights_method.lower() if wm == "queen": w = libpysal.weights.Queen.from_dataframe(gdf, use_index=True) elif wm == "rook": w = libpysal.weights.Rook.from_dataframe(gdf, use_index=True) elif wm == "distance": thr = distance_threshold if target_crs.upper() == "EPSG:4326": thr = distance_threshold / 111000.0 # meters -> degrees w = libpysal.weights.DistanceBand.from_dataframe( gdf, threshold=thr, binary=False, use_index=True ) else: return {"status":"error","message":f"Unknown weights_method: {weights_method}"} w.transform = "r" # handle islands by dropping rows and rebuilding weights if w.islands: keep_idx = [i for i in range(gdf.shape[0]) if i not in set(w.islands)] if not keep_idx: return {"status":"error","message":"All units are islands under current weights; adjust weights_method/threshold."} gdf = gdf.iloc[keep_idx].reset_index(drop=True) Y = Y[keep_idx, :] if wm == "queen": w = libpysal.weights.Queen.from_dataframe(gdf, use_index=True) elif wm == "rook": w = libpysal.weights.Rook.from_dataframe(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(gdf, threshold=thr, binary=False, use_index=True) w.transform = "r" # --- Spatial Markov --- try: from giddy.markov import Spatial_Markov except ModuleNotFoundError: return {"status": "error", "message": "The 'giddy' package is not installed. Install with: pip install giddy"} sm = Spatial_Markov( Y, w, k=int(k), m=int(m), permutations=int(permutations), fixed=bool(fixed), discrete=False, fill_empty_classes=bool(fill_empty_classes), variable_name=";".join(value_cols) ) # --- package results (JSON-safe) --- def tolist(x): """Recursively convert numpy arrays and nested structures to lists.""" if x is None: return None try: if isinstance(x, np.ndarray): return x.tolist() elif isinstance(x, (list, tuple)): return [tolist(xx) for xx in x] elif isinstance(x, dict): return {k: tolist(v) for k, v in x.items()} elif isinstance(x, (np.integer, np.floating)): return float(x) if isinstance(x, np.floating) else int(x) else: return x except Exception: # Fallback: try to convert to native Python type try: return float(x) if isinstance(x, (np.floating, float)) else int(x) if isinstance(x, (np.integer, int)) else str(x) except Exception: return x # Safer preview: keep geometry, write WKT to a new column, then drop geometry preview = gdf[[*value_cols, "geometry"]].head(5).copy() preview["geometry_wkt"] = preview["geometry"].apply( lambda g: g.wkt if g is not None and not g.is_empty else None ) preview = pd.DataFrame(preview.drop(columns="geometry")).to_dict(orient="records") result = { "n_regions": int(Y.shape[0]), "n_periods": int(Y.shape[1]), "k_classes_y": int(k), "m_classes_lag": int(m), "weights_method": weights_method.lower(), "value_columns": value_cols, "discretization": { "cutoffs_y": tolist(getattr(sm, "cutoffs", None)), "cutoffs_lag": tolist(getattr(sm, "lag_cutoffs", None)), "fixed": bool(fixed) }, "global_transition_prob_p": tolist(sm.p), # (k x k) "conditional_transition_prob_P": tolist(sm.P), # (m x k x k) "global_steady_state_s": tolist(sm.s), # (k,) "conditional_steady_states_S": tolist(sm.S), # (m x k) "tests": { "chi2_total_x2": float(getattr(sm, "x2", np.nan)), "chi2_df": int(getattr(sm, "x2_dof", -1)), "chi2_pvalue": float(getattr(sm, "x2_pvalue", np.nan)), "Q": float(getattr(sm, "Q", np.nan)), "Q_p_value": float(getattr(sm, "Q_p_value", np.nan)), "LR": float(getattr(sm, "LR", np.nan)), "LR_p_value": float(getattr(sm, "LR_p_value", np.nan)), }, "data_preview": preview, } msg = "Spatial Markov completed successfully" if wm == "distance" and target_crs.upper() == "EPSG:4326": msg += f" (threshold interpreted as {distance_threshold/111000.0:.6f} degrees)." # Apply tolist conversion recursively to the entire result result = tolist(result) return {"status": "success", "message": msg, "result": result} except Exception as e: return {"status": "error", "message": f"Failed to run Spatial Markov: {str(e)}"}