Add batch PLSS queries with optional county validation

- get_plss_details_batch: query multiple coordinates in parallel
- Census TIGERweb integration for county lookups (layer 82)
- Returns unique_townships for discovering T/R combinations in a grid
- County data useful for validating sample points fall within expected area

Version bump to 2026.01.25

README.md

@@ -52,10 +52,12 @@ uvx mcblmplss
 claude mcp add blm "uvx mcblmplss"
 ```
 
-## Tools (18 total)
+## Tools (19 total)
 
 Each data type has two tools: a human-readable version and a `_details` version returning structured data.
 
+Plus a batch tool for efficient multi-coordinate queries.
+
 ### PLSS — Public Land Survey System
 
 ```
@@ -66,6 +68,38 @@ State: CO
 PLSS ID: CO060010S0680W0SN090
 ```
 
+### Batch PLSS Queries
+
+Query multiple coordinates in parallel with optional county validation:
+
+```
+> get_plss_details_batch(
+    coordinates=[[44.5, -115.8], [44.6, -115.9], [44.7, -116.0]],
+    include_county=True
+  )
+
+{
+  "results": [
+    {
+      "latitude": 44.5,
+      "longitude": -115.8,
+      "township": {"township": "13N", "range": "5E", ...},
+      "section": {"section": "30", ...},
+      "county": {"name": "Valley County", "fips": "16085", ...}
+    },
+    ...
+  ],
+  "total": 3,
+  "successful": 3,
+  "failed": 0,
+  "unique_townships": ["13N 5E", "14N 4E"]
+}
+```
+
+The `unique_townships` field is useful for discovering which Township/Range combinations exist within a sampling grid.
+
+**County data** comes from Census TIGERweb (layer 82), not BLM PLSS. Use `include_county=True` when you need to validate that sampled points fall within an expected county.
+
 ### Land Manager — Surface Management Agency
 
 ```
@@ -177,11 +211,12 @@ Data is available for **30 states** where the Public Land Survey System was used
 
 ## Data Sources
 
-All data comes from official BLM ArcGIS REST services:
+All data comes from official government ArcGIS REST services (BLM + Census TIGERweb):
 
 | Data | Source | Typical Update |
 |------|--------|----------------|
 | PLSS | [BLM National PLSS CadNSDI](https://gis.blm.gov/arcgis/rest/services/Cadastral/BLM_Natl_PLSS_CadNSDI/MapServer) | Quarterly |
+| County | [Census TIGERweb](https://tigerweb.geo.census.gov/arcgis/rest/services/TIGERweb/tigerWMS_Current/MapServer) | Annual |
 | Surface Management | [BLM SMA](https://gis.blm.gov/arcgis/rest/services/lands/BLM_Natl_SMA_LimitedScale/MapServer) | Annual |
 | Mining Claims | [BLM MLRS](https://gis.blm.gov/nlsdb/rest/services/Mining_Claims/MiningClaims/MapServer) | Weekly |
 | Grazing | [BLM Grazing Allotment](https://gis.blm.gov/arcgis/rest/services/range/BLM_Natl_Grazing_Allotment/MapServer) | Annual |

pyproject.toml

@@ -1,6 +1,6 @@
 [project]
 name = "mcblmplss"
-version = "2025.01.25"
+version = "2026.01.25"
 description = "MCP server for querying BLM land data: PLSS, land manager, mining claims, grazing, wild horses, recreation, wilderness, rivers, and ACECs"
 readme = "README.md"
 license = "MIT"

src/mcblmplss/census.py

@@ -0,0 +1,120 @@
+"""
+Census Bureau TIGERweb client for county lookups.
+
+Provides county information via the Census Bureau's TIGERweb ArcGIS service.
+This complements BLM data which doesn't include county boundaries.
+"""
+
+import httpx
+from pydantic import BaseModel, Field
+
+# Census TIGERweb MapServer
+CENSUS_URL = "https://tigerweb.geo.census.gov/arcgis/rest/services/TIGERweb/tigerWMS_Current/MapServer"
+LAYER_COUNTY = 82
+
+
+class CountyInfo(BaseModel):
+    """County information from Census TIGERweb."""
+
+    name: str = Field(..., description="County name (e.g., 'Valley County')")
+    fips: str = Field(..., description="Full FIPS code: state + county (e.g., '16085')")
+    county_fips: str = Field(..., description="County FIPS code only (e.g., '085')")
+    state_fips: str = Field(..., description="State FIPS code (e.g., '16' for Idaho)")
+
+
+class CensusAPIError(Exception):
+    """Error communicating with Census TIGERweb services."""
+
+    pass
+
+
+async def get_county(latitude: float, longitude: float, timeout: float = 30.0) -> CountyInfo | None:
+    """
+    Get county information for coordinates.
+
+    Args:
+        latitude: WGS84 latitude
+        longitude: WGS84 longitude
+        timeout: Request timeout in seconds
+
+    Returns:
+        CountyInfo if found, None if no county data available
+
+    Raises:
+        CensusAPIError: If the API request fails
+    """
+    params = {
+        "f": "json",
+        "geometry": f"{longitude},{latitude}",
+        "geometryType": "esriGeometryPoint",
+        "sr": "4326",
+        "layers": f"all:{LAYER_COUNTY}",
+        "tolerance": "1",
+        "mapExtent": f"{longitude - 1},{latitude - 1},{longitude + 1},{latitude + 1}",
+        "imageDisplay": "100,100,96",
+        "returnGeometry": "false",
+    }
+
+    try:
+        async with httpx.AsyncClient(timeout=timeout) as client:
+            response = await client.get(f"{CENSUS_URL}/identify", params=params)
+            response.raise_for_status()
+            data = response.json()
+    except httpx.TimeoutException:
+        raise CensusAPIError("Census API request timed out.")
+    except httpx.HTTPStatusError as e:
+        raise CensusAPIError(f"Census API returned error {e.response.status_code}")
+    except httpx.RequestError as e:
+        raise CensusAPIError(f"Failed to connect to Census API: {e}")
+
+    if "error" in data:
+        error_msg = data["error"].get("message", "Unknown error")
+        raise CensusAPIError(f"Census API error: {error_msg}")
+
+    results = data.get("results", [])
+    if not results:
+        return None
+
+    # Find the county layer result
+    for result in results:
+        if result.get("layerId") == LAYER_COUNTY:
+            attrs = result.get("attributes", {})
+            name = attrs.get("NAME", "")
+            geoid = attrs.get("GEOID", "")
+            county_fips = attrs.get("COUNTY", "")
+
+            if geoid and len(geoid) >= 5:
+                return CountyInfo(
+                    name=name,
+                    fips=geoid,
+                    county_fips=county_fips,
+                    state_fips=geoid[:2],
+                )
+
+    return None
+
+
+async def get_counties_batch(
+    coordinates: list[tuple[float, float]],
+    timeout: float = 30.0,
+) -> list[CountyInfo | None]:
+    """
+    Get county information for multiple coordinates in parallel.
+
+    Args:
+        coordinates: List of (latitude, longitude) tuples
+        timeout: Request timeout per coordinate
+
+    Returns:
+        List of CountyInfo or None for each coordinate
+    """
+    import asyncio
+
+    async def safe_get_county(lat: float, lon: float) -> CountyInfo | None:
+        try:
+            return await get_county(lat, lon, timeout)
+        except CensusAPIError:
+            return None
+
+    tasks = [safe_get_county(lat, lon) for lat, lon in coordinates]
+    return await asyncio.gather(*tasks)

src/mcblmplss/client.py

@@ -4,6 +4,7 @@ Shared HTTP client for BLM ArcGIS REST API queries.
 Provides common identify/query operations against BLM MapServer endpoints.
 """
 
+import asyncio
 from typing import Any
 
 import httpx
@@ -143,6 +144,34 @@ class BLMClient:
         features = data.get("features", [])
         return [f.get("attributes", {}) for f in features]
 
+    async def identify_batch(
+        self,
+        base_url: str,
+        coordinates: list[tuple[float, float]],
+        layers: str,
+        tolerance: int = 1,
+        return_geometry: bool = False,
+    ) -> list[list[dict[str, Any]] | BLMAPIError]:
+        """
+        Perform identify operations on multiple coordinates in parallel.
+
+        Args:
+            base_url: MapServer URL (e.g., https://gis.blm.gov/.../MapServer)
+            coordinates: List of (latitude, longitude) tuples
+            layers: Layer specification (e.g., "all:1,2" or "visible:0")
+            tolerance: Pixel tolerance for hit detection
+            return_geometry: Whether to include geometry in response
+
+        Returns:
+            List of results for each coordinate. Each result is either a list
+            of result dictionaries, or a BLMAPIError if that request failed.
+        """
+        tasks = [
+            self.identify(base_url, lat, lon, layers, tolerance, return_geometry)
+            for lat, lon in coordinates
+        ]
+        return await asyncio.gather(*tasks, return_exceptions=True)
+
 
 # Shared client instance
 blm_client = BLMClient(timeout=60.0)

src/mcblmplss/mixins/__init__.py

@@ -4,10 +4,11 @@ MCP Mixins for BLM data services.
 Each mixin provides tools for a specific BLM data domain.
 """
 
+from mcblmplss.census import CountyInfo
 from mcblmplss.mixins.acec import ACECMixin
 from mcblmplss.mixins.grazing import GrazingMixin
 from mcblmplss.mixins.mining_claims import MiningClaimsMixin
-from mcblmplss.mixins.plss import PLSSMixin
+from mcblmplss.mixins.plss import BatchPLSSResult, PLSSMixin
 from mcblmplss.mixins.recreation import RecreationMixin
 from mcblmplss.mixins.surface_management import SurfaceManagementMixin
 from mcblmplss.mixins.wild_horses import WildHorseMixin
@@ -15,6 +16,7 @@ from mcblmplss.mixins.wild_rivers import WildRiversMixin
 from mcblmplss.mixins.wilderness import WildernessMixin
 
 __all__ = [
+    # Mixins
     "ACECMixin",
     "PLSSMixin",
     "SurfaceManagementMixin",
@@ -24,4 +26,7 @@ __all__ = [
     "WildRiversMixin",
     "WildernessMixin",
     "RecreationMixin",
+    # Models
+    "BatchPLSSResult",
+    "CountyInfo",
 ]

src/mcblmplss/mixins/plss.py

@@ -7,6 +7,7 @@ Provides tools for querying Section, Township, and Range from coordinates.
 from fastmcp.contrib.mcp_mixin import MCPMixin, mcp_tool
 from pydantic import BaseModel, Field
 
+from mcblmplss.census import CountyInfo, get_counties_batch
 from mcblmplss.client import BLMAPIError, blm_client
 
 # BLM Cadastral MapServer
@@ -34,9 +35,26 @@ class PLSSResult(BaseModel):
     longitude: float
     township: PLSSLocation | None = None
     section: PLSSLocation | None = None
+    county: CountyInfo | None = Field(
+        None, description="County information from Census TIGERweb"
+    )
     error: str | None = None
 
 
+class BatchPLSSResult(BaseModel):
+    """Result from batch PLSS coordinate query."""
+
+    results: list[PLSSResult] = Field(
+        ..., description="Individual results for each coordinate"
+    )
+    total: int = Field(..., description="Total number of coordinates queried")
+    successful: int = Field(..., description="Number of successful lookups")
+    failed: int = Field(..., description="Number of failed lookups")
+    unique_townships: list[str] = Field(
+        ..., description="Unique Township/Range combinations found (e.g., '13N 5E')"
+    )
+
+
 def _parse_township(attrs: dict) -> PLSSLocation:
     """Parse township attributes from API response."""
     twp_num = attrs.get("TWNSHPNO", "").lstrip("0") or "0"
@@ -179,3 +197,119 @@ class PLSSMixin(MCPMixin):
             township=township_data,
             section=section_data,
         )
+
+    @mcp_tool(
+        name="get_plss_details_batch",
+        description=(
+            "Get PLSS details for multiple coordinates in parallel. "
+            "Returns structured data with unique township/range summary."
+        ),
+    )
+    async def get_plss_details_batch(
+        self,
+        coordinates: list[tuple[float, float]] = Field(
+            description="List of [latitude, longitude] coordinate pairs"
+        ),
+        include_county: bool = Field(
+            default=False,
+            description="Include county lookup from Census TIGERweb (validates location)",
+        ),
+    ) -> BatchPLSSResult:
+        """
+        Get PLSS details for multiple coordinates in parallel.
+
+        Optimized for batch queries (e.g., grid sampling a county). Returns
+        results with unique township/range summary for T/R discovery workflows.
+
+        Args:
+            coordinates: List of (latitude, longitude) tuples
+            include_county: If True, also query Census TIGERweb for county info
+
+        Returns:
+            BatchPLSSResult with individual results and unique T/R summary
+        """
+        # Query PLSS data in parallel
+        batch_results = await blm_client.identify_batch(
+            PLSS_URL,
+            coordinates,
+            f"all:{LAYER_TOWNSHIP},{LAYER_SECTION}",
+        )
+
+        # Optionally query county data in parallel
+        counties: list[CountyInfo | None] = []
+        if include_county:
+            counties = await get_counties_batch(coordinates)
+        else:
+            counties = [None] * len(coordinates)
+
+        # Process results
+        results: list[PLSSResult] = []
+        unique_tr: set[str] = set()
+        successful = 0
+        failed = 0
+
+        for i, (lat, lon) in enumerate(coordinates):
+            api_result = batch_results[i]
+            county_info = counties[i] if include_county else None
+
+            # Handle API errors
+            if isinstance(api_result, Exception):
+                results.append(
+                    PLSSResult(
+                        latitude=lat,
+                        longitude=lon,
+                        county=county_info,
+                        error=str(api_result),
+                    )
+                )
+                failed += 1
+                continue
+
+            # Parse PLSS data
+            if not api_result:
+                results.append(
+                    PLSSResult(
+                        latitude=lat,
+                        longitude=lon,
+                        county=county_info,
+                        error="No PLSS data found. Location may be outside surveyed areas.",
+                    )
+                )
+                failed += 1
+                continue
+
+            township_data: PLSSLocation | None = None
+            section_data: PLSSLocation | None = None
+
+            for result in api_result:
+                layer_id = result.get("layerId")
+                attrs = result.get("attributes", {})
+
+                if layer_id == LAYER_TOWNSHIP and township_data is None:
+                    township_data = _parse_township(attrs)
+                elif layer_id == LAYER_SECTION and section_data is None:
+                    section_data = _parse_section(attrs, township_data)
+
+            # Track unique T/R combinations
+            if township_data:
+                tr_key = f"{township_data.township} {township_data.range}"
+                unique_tr.add(tr_key)
+
+            results.append(
+                PLSSResult(
+                    latitude=lat,
+                    longitude=lon,
+                    township=township_data,
+                    section=section_data,
+                    county=county_info,
+                )
+            )
+            successful += 1
+
+        return BatchPLSSResult(
+            results=results,
+            total=len(coordinates),
+            successful=successful,
+            failed=failed,
+            unique_townships=sorted(unique_tr),
+        )