GIS MCP Server

Instructions

Implementation Reference

• src/gis_mcp/pysal_functions.py:1322-1461 (handler)

  The handler function for 'dynamic_lisa' tool, implementing dynamic Local Indicators of Spatial Association using the giddy.directional.Rose class. Computes directional statistics between two time periods in shapefile data.
  @gis_mcp.tool() def dynamic_lisa( shapefile_path: str, value_columns: Union[str, List[str]], # exactly two columns: [t0, t1] target_crs: str = "EPSG:4326", weights_method: str = "queen", # 'queen'|'rook'|'distance' distance_threshold: float = 100000, # meters (converted to degrees if 4326) k: int = 8, # number of rose sectors permutations: int = 99, # 0 = skip inference alternative: str = "two.sided", # 'two.sided'|'positive'|'negative' relative: bool = True, # divide each column by its mean drop_na: bool = True # drop features with NA in either column ) -> Dict[str, Any]: """Run dynamic LISA (directional LISA) with giddy.directional.Rose. Returns sector counts, angles, vector lengths, and (optionally) permutation p-values. """ try: # --- sanitize --- for s in ("shapefile_path","target_crs","weights_method","alternative"): v = locals()[s] if isinstance(v, str): locals()[s] = v.replace("`","") if isinstance(value_columns, str): cols = [c.strip() for c in value_columns.replace("`","").split(",") if c.strip()] else: cols = list(value_columns) if not os.path.exists(shapefile_path): return {"status":"error","message":f"Shapefile not found: {shapefile_path}"} if len(cols) != 2: return {"status":"error","message":"value_columns must be exactly two columns: [start_time, end_time]."} # --- load + project --- gdf = gpd.read_file(shapefile_path) missing = [c for c in cols if c not in gdf.columns] if missing: return {"status":"error","message":f"Columns not found: {missing}"} gdf = gdf.to_crs(target_crs) # --- prepare Y (n x 2) --- Y = gdf[cols].astype(float).to_numpy(copy=True) if drop_na: mask = ~np.any(np.isnan(Y), axis=1) gdf = gdf.loc[mask].reset_index(drop=True) Y = Y[mask, :] 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) else: return {"status":"error","message":f"Unknown weights_method: {weights_method}"} w.transform = "r" # drop islands (units with no neighbors) if w.islands: keep = [i for i in range(gdf.shape[0]) if i not in set(w.islands)] if not keep: return {"status":"error","message":"All units are islands under current weights."} gdf = gdf.iloc[keep].reset_index(drop=True) Y = Y[keep, :] 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) w.transform = "r" # --- Dynamic LISA (Rose) --- from giddy.directional import Rose rose = Rose(Y, w, k=int(k)) # computes cuts, counts, theta, r, bins, lag internally # permutation inference (optional) pvals = None expected_perm = larger_perm = smaller_perm = None if permutations and permutations > 0: rose.permute(permutations=int(permutations), alternative=alternative) pvals = np.asarray(getattr(rose, "p", None)).tolist() if hasattr(rose, "p") else None expected_perm = np.asarray(getattr(rose, "expected_perm", [])).tolist() larger_perm = np.asarray(getattr(rose, "larger_perm", [])).tolist() smaller_perm = np.asarray(getattr(rose, "smaller_perm", [])).tolist() # preview def _wkt_safe(g): try: return g.wkt except Exception: return str(g) preview = gdf[[*cols, gdf.geometry.name]].head(5).copy() preview[gdf.geometry.name] = preview[gdf.geometry.name].apply(_wkt_safe) preview = preview.rename(columns={gdf.geometry.name: "geometry_wkt"}).to_dict(orient="records") result = { "n_regions": int(Y.shape[0]), "k_sectors": int(k), "weights_method": wm, "value_columns": cols, "cuts_radians": np.asarray(rose.cuts).tolist(), # sector edges "sector_counts": np.asarray(rose.counts).tolist(), # count per sector "angles_theta_rad": np.asarray(rose.theta).tolist(), # one per region "vector_lengths_r": np.asarray(rose.r).tolist(), # one per region "bins_used": np.asarray(rose.bins).tolist(), "inference": { "permutations": int(permutations), "alternative": alternative, "p_values_by_sector": pvals, "expected_counts_perm": expected_perm, "larger_or_equal_counts": larger_perm, "smaller_or_equal_counts": smaller_perm }, "data_preview": preview } msg = "Dynamic LISA (Rose) completed successfully" if wm == "distance" and target_crs.upper() == "EPSG:4326": msg += f" (threshold interpreted as {distance_threshold/111000.0:.6f} degrees)." return {"status":"success","message":msg,"result":result} except Exception as e: return {"status":"error","message":f"Failed to run Dynamic LISA: {str(e)}"}
• src/gis_mcp/__init__.py:9-10 (registration)

  Imports pysal_functions.py, which triggers registration of the dynamic_lisa tool via @gis_mcp.tool() decorator.
  from .pysal_functions import * from .save_tool import *