GIS MCP Server

Instructions

Implementation Reference

• src/gis_mcp/pysal_functions.py:396-453 (handler)

  Main handler function for the 'join_counts_local' tool. Loads spatial data, creates distance-based weights, handles isolated units, computes local join counts using esda.Join_Counts_Local, and returns results with preview.

  @gis_mcp.tool()
  def join_counts_local(shapefile_path: str, dependent_var: str = "LAND_USE", target_crs: str = "EPSG:4326",
                        distance_threshold: float = 100000) -> Dict[str, Any]:
      """Local Join Counts."""
      gdf, y, w, (threshold, unit), err = pysal_load_data(shapefile_path, dependent_var, target_crs, distance_threshold)
      if err:
          return {"status": "error", "message": err}
  
      # Handle islands - if all points are islands, fall back to KNN weights for connectivity
      import libpysal
      if w.islands:
          if len(w.islands) == len(gdf):
              # All points are islands - fall back to KNN weights
              try:
                  # Use k=4 for a 5x5 grid to ensure connectivity
                  w = libpysal.weights.KNN.from_dataframe(gdf, k=4)
                  w.transform = 'r'
              except Exception as e:
                  return {"status": "error", "message": f"All units are islands and KNN fallback failed: {str(e)}"}
          else:
              # Some islands - filter them out
              keep_idx = [i for i in range(len(gdf)) if i not in set(w.islands)]
              if len(keep_idx) == 0:
                  return {"status": "error", "message": "All units are islands (no neighbors). Try increasing distance_threshold."}
              # Filter data
              gdf_filtered = gdf.iloc[keep_idx].reset_index(drop=True)
              y_filtered = y[keep_idx]
              # Rebuild weights without islands using the same threshold
              w_filtered = libpysal.weights.DistanceBand.from_dataframe(
                  gdf_filtered, 
                  threshold=threshold,  # Use the effective threshold already calculated in pysal_load_data
                  binary=False
              )
              w_filtered.transform = 'r'
              gdf, y, w = gdf_filtered, y_filtered, w_filtered
  
      import esda
      stat = esda.Join_Counts_Local(y, w)
      preview = gdf[['geometry', dependent_var]].head(5).copy()
      preview['geometry'] = preview['geometry'].apply(lambda g: g.wkt)
  
      # Join_Counts_Local has LJC attribute
      ljc_val = None
      if hasattr(stat, "LJC"):
          ljc_val = stat.LJC.tolist() if hasattr(stat.LJC, "tolist") else list(stat.LJC)
      elif hasattr(stat, "local_join_counts"):
          ljc_val = stat.local_join_counts.tolist() if hasattr(stat.local_join_counts, "tolist") else list(stat.local_join_counts)
      elif hasattr(stat, "ljc"):
          ljc_val = stat.ljc.tolist() if hasattr(stat.ljc, "tolist") else list(stat.ljc)
      
      return {
          "status": "success",
          "message": f"Local Join Counts completed successfully (threshold: {threshold} {unit})",
          "result": {
              "local_join_counts": ljc_val,
              "data_preview": preview.to_dict(orient="records")
          }
      }

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

  Helper function used by join_counts_local (and other tools) to load GeoDataFrame from shapefile, reproject, extract dependent variable, create and row-standardize distance band spatial weights, and handle basic island isolation.

  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

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

  GIS resource endpoint listing available ESDA/PySAL operations, including 'join_counts_local' as one of the supported tools.

  @gis_mcp.resource("gis://operations/esda")
  def get_spatial_operations() -> Dict[str, List[str]]:
      """List available spatial analysis operations. This is for esda library. They are using pysal library."""
      return {
          "operations": [
              "getis_ord_g",
              "morans_i",
              "gearys_c",
              "gamma_statistic",
              "moran_local",
              "getis_ord_g_local",
              "join_counts",
              "join_counts_local",
              "adbscan"
          ]
      }