georest.esri.edw

USFS Enterprise Data Warehouse (EDW) REST API client.

Provides search, metadata inspection, and spatial feature queries against the ArcGIS REST services at https://apps.fs.usda.gov/arcx/rest/services/EDW.

Quick start:

from georest.esri import edw

# Search for fire-related services (EDW catalog only - see search_edw_services)
services = edw.search_edw_services("fire")

# Get layer info
info = edw.get_service_info("EDW_MTBS_01")

# Query features as GeoJSON (layer 63 = "Burned Area Boundaries (All Years)";
# the lower-numbered layers are single years)
geojson = edw.query_features("EDW_MTBS_01", 63,
    where="fire_name LIKE '%CAMERON PEAK%'",
    out_fields="fire_name,acres,year")

Field naming gotcha: where and out_fields take field names, which on most EDW layers are lowercase (fire_name, acres, year). The uppercase forms shown in ArcGIS clients (FIRE_NAME, ACRES, YEAR) are display aliases; passing those in out_fields makes the server reject the whole query with “Failed to execute query.” Use get_layer_info(service, layer)[“fields”] to get the real names.

Copyright 2026 Ryan Rock and Ian Housman

Licensed under the Apache License, Version 2.0 (the “License”); you may not use this file except in compliance with the License. You may obtain a copy of the License at

http://www.apache.org/licenses/LICENSE-2.0

Module Attributes

VALID_THEMES	Themes considered valid for a _SERVICE_THEMES entry, and the set search_edw_services(theme=...) accepts.
LAYER_ROLES	Roles assigned by get_layer_roles().

Functions

get_detail_layers(service_name)	An EDW service's queryable layers, minus coarse cartographic duplicates.
get_layer_info(service_name, layer_id)	Get detailed metadata for a specific layer: fields, geometry type, capabilities.
get_layer_metadata(service_name, layer_id[, ...])	Fetch and parse a layer's FGDC/ISO metadata XML (the /metadata operation).
get_layer_roles(service_name)	Every layer in an EDW service, tagged with its cartographic role.
get_service_info(service_name)	Get metadata for an EDW MapServer service: description, layers, spatial ref.
query_features(service_name, layer_id[, ...])	Query features from an EDW layer, optionally filtered by spatial intersection.
query_features_analytic(service_name, ...[, ...])	Run SQL window-function analytics (rank, running total, LAG/LEAD, ...) against an EDW layer via the ArcGIS queryAnalytic operation.
query_features_with_pagination(service_name, ...)	Query features with automatic pagination past the server's per-request cap.
search_edw_services([query, theme])	Search the USFS EDW catalog ONLY by keyword and/or theme.
theme_drift(services)	Compare the live EDW catalog against the hand-maintained theme table.
top_n_per_group(service_name, layer_id, ...)	Get the top N features by field within each group_by partition.

georest.esri.edw.VALID_THEMES = frozenset({'biota', 'boundaries', 'environment', 'geoscientific', 'inland_waters', 'planning_cadastre', 'structure', 'transportation'})

Themes considered valid for a _SERVICE_THEMES entry, and the set search_edw_services(theme=…) accepts.

georest.esri.edw.theme_drift(services) → dict[str, list[str]]

Compare the live EDW catalog against the hand-maintained theme table.

_SERVICE_THEMES powers theme filtering and description matching in search_edw_services, and is curated by hand - so it silently rots as EDW publishes and retires services.

Only MapServer services are expected to carry a theme. The catalog also publishes the odd GPServer (a geoprocessing tool, not a queryable data service - currently RAVG_DataExtract_01), which has no layers and nothing to categorise; those are reported under “non_mapserver” for visibility rather than counted as drift.

Takes the list of dicts from search_edw_services. Returns:

untyped - live MapServer services with no theme entry. They still

appear in search results, but as theme “uncategorized” and with no description to match keywords against.

orphaned - theme entries whose service is gone from the catalog. bad_theme / missing_desc - malformed entries. non_mapserver - informational; not drift.

Example:

from georest.esri import edw

# Compare the live catalog with georest's theme table (what
# georest-check-themes reports)
edw.theme_drift(edw.search_edw_services())

georest.esri.edw.search_edw_services(query: str = '', theme: str = '') → list[dict[str, str]]

Search the USFS EDW catalog ONLY by keyword and/or theme.

IMPORTANT - scope: this searches exactly one flat catalog, EDW_BASE_URL (https://apps.fs.usda.gov/arcx/rest/services/EDW). It does NOT search IIPP, ArcGIS Online, or any other portal. To search IIPP (or another portal) instead, use georest.esri.portal.searchPortal(query, portal="iipp") - that function hits the portal’s /sharing/rest/search endpoint directly and has no EDW-specific theme/keyword-alias matching.

Uses three matching strategies (results are deduplicated): 1. Substring match on service name (e.g. “mtbs” → EDW_MTBS_01) 2. Keyword alias expansion (e.g. “riparian” → inland_waters theme + stream services) 3. Theme filter (e.g. theme=”inland_waters”)

Parameters:

• query – Search keyword (case-insensitive). Matches against service names, theme descriptions, and keyword aliases.

• theme – Filter by theme category. Valid themes: biota, boundaries, environment, geoscientific, inland_waters, planning_cadastre, structure, transportation. Pass “” to skip theme filtering.

Returns a list of dicts with keys: name, type, url, theme, description. If both query and theme are empty, returns all services.

Example:

from georest.esri import edw

for s in edw.search_edw_services("fire"):          # by keyword
    print(s["name"], "-", s["description"])

edw.search_edw_services(theme="recreation")         # or by theme

georest.esri.edw.get_service_info(service_name: str) → dict[str, Any]

Get metadata for an EDW MapServer service: description, layers, spatial ref.

Parameters:

service_name – Short service name, e.g. “EDW_MTBS_01”.

Returns:

name, description, spatialReference, fullExtent, layers (id, name, defaultVisibility, minScale, maxScale).

Note there is no geometryType here: the MapServer catalog’s layers array doesn’t carry it. Use get_layer_info(service_name, layer_id) for a layer’s geometryType and fields.

Return type:

Dict with keys

Example:

from georest.esri import edw

info = edw.get_service_info("EDW_MTBS_01")
[(layer["id"], layer["name"]) for layer in info["layers"]]

georest.esri.edw.LAYER_ROLES = ('group', 'detail', 'coarse', 'standalone')

Roles assigned by get_layer_roles(). See that function for the criteria.

georest.esri.edw.get_layer_roles(service_name: str) → list[dict[str, Any]]

Every layer in an EDW service, tagged with its cartographic role.

~18% of EDW services (19 of the 105 whose layer trees were surveyed - EDW_NorWeST_StreamTemperatures_01, the EDW_HydroFlowMetrics* family, EDW_County_01, EDW_State_01, EDW_RangerDistricts_01/03, …) publish a theme more than once: a coarse layer meant only for small-scale cartographic overview (drawn when zoomed far out, switched off once you zoom in) alongside a detailed layer carrying the real geometry/attribute resolution (drawn only once zoomed in). For feature-level work (query_features, spatial joins, attribute stats) the coarse layer is almost never what you want - querying it can silently return generalized geometry rather than an error.

Rather than drop anything, this returns the full layer list with each entry tagged, so a mis-grouping is visible instead of silently losing a layer. Each dict is get_service_info()[‘layers’] shape plus:

type, parentLayerId - as reported by the server role - one of:

“group” a Group Layer container. Has no /query endpoint,

so it is never directly queryable.

“detail” the full-resolution member of a coarse/detail split “coarse” a small-scale sibling superseded by a “detail”

layer in the same group

“standalone” no split applies to its group; use it as-is

groupKey - the (“parent”, id) or (“name”, base) key its siblings

were matched on (None for group layers), so you can audit why two layers were treated as versions of each other

Role criteria, within one group: a member is coarse if maxScale > 0 (i.e. it stops drawing past some scale), and detail if minScale > 0 and maxScale == 0. A group only splits when it has at least one of each; otherwise every member is “standalone”. Note a coarse layer may itself have minScale > 0 (EDW_CongressionalDistricts_04 publishes its national tier as minScale=77000000, maxScale=1000000), which is why the coarse test is maxScale > 0 alone. Any member of a splitting group that is not detail - including a mid-resolution tier banded on both sides - is tagged coarse.

Grouping caveat: siblings are matched by shared parent Group Layer when there is one, else by name with a trailing “- National/Regional Extent” or “(National/Regional Extent)” stripped. Parent-based matching assumes a group holds one dataset at several resolutions. That holds throughout the current EDW catalog, but a group mixing genuinely different datasets could see one tagged “coarse” against another - inspect groupKey if a result looks wrong.

Raises:

RuntimeError – if the service errors, or returns no layers at all.

Example:

from georest.esri import edw

for layer in edw.get_layer_roles("EDW_ForestSystemBoundaries_01"):
    print(layer["id"], layer["role"], layer["name"])

georest.esri.edw.get_detail_layers(service_name: str) → list[dict[str, Any]]

An EDW service’s queryable layers, minus coarse cartographic duplicates.

Convenience filter over get_layer_roles(): keeps the “detail” and “standalone” layers and sets aside Group Layers (not queryable) and “coarse” small-scale duplicates. Safe to call on any service - one with no coarse/detail split returns all of its data layers.

Use get_layer_roles() directly when you need to see what was set aside and why; this function reports only what survived.

Example:

from georest.esri import edw

[layer["name"] for layer in edw.get_detail_layers("EDW_ForestSystemBoundaries_01")]

georest.esri.edw.get_layer_info(service_name: str, layer_id: int) → dict[str, Any]

Get detailed metadata for a specific layer: fields, geometry type, capabilities.

Parameters:

• service_name – Short service name, e.g. “EDW_MTBS_01”.

• layer_id – Layer ID within the service.

Returns:

name, geometryType, description, fields, extent, maxRecordCount, supportedQueryFormats, capabilities (list of supported operations, e.g. [“Map”, “Query”, “Data”]), advancedQueryCapabilities (dict of supportsX booleans, e.g. supportsQueryAnalytic, supportsStatistics - varies by service).

Return type:

Dict with keys

Example:

from georest.esri import edw

layer = edw.get_layer_info("EDW_MTBS_01", 63)   # Burned Area Boundaries (All Years)
[f["name"] for f in layer["fields"]]

georest.esri.edw.get_layer_metadata(service_name: str, layer_id: int, include_domains: bool = False) → dict[str, Any]

Fetch and parse a layer’s FGDC/ISO metadata XML (the /metadata operation).

This is a separate document from get_layer_info’s REST JSON - it carries things the layer JSON doesn’t: a dataset abstract/purpose, keywords, and per-field definitions (e.g. what “Post_ID” actually means).

Parameters:

• service_name – Short service name, e.g. “EDW_MTBS_01”.

• layer_id – Layer ID within the service.

• include_domains – If True, add a domain key to each attribute with its enumerated/range domain (see attributes below). Free-text (“udom”) domains are not structured data and are ignored.

Returns:

title, abstract, purpose, credit (str, “” if not present) keywords (sorted list of theme keywords) attributes: list of {name, alias, type, definition,

definition_source} per field. definition is “” for fields the data provider didn’t document (e.g. OBJECTID clones, raw thresholds). If include_domains is True, each entry also has domain: None, or a dict with either/both:

• ”coded_values”: [{“value”, “definition”, “definition_source”}, …] for enumerated domains

• ”range”: {“min”, “max”} for numeric range domains

A field can have both (e.g. a numeric field using -9999 as a documented “no data” sentinel alongside its real range).

Return type:

Dict with keys

Raises:

RuntimeError – If the request fails or the response isn’t valid XML (some layers have no metadata document configured).

Example:

from georest.esri import edw

meta = edw.get_layer_metadata("EDW_MTBS_01", 63)
meta["title"], meta["abstract"]

georest.esri.edw.query_features(service_name: str, layer_id: int, geometry: dict | str | None = None, geometry_type: str = 'esriGeometryEnvelope', spatial_rel: str = 'esriSpatialRelIntersects', where: str = '1=1', out_fields: str = '*', max_features: int = 1000, out_sr: int = 4326, return_count_only: bool = False, in_sr: int | None = None) → dict

Query features from an EDW layer, optionally filtered by spatial intersection.

Parameters:

• service_name – Short service name, e.g. “EDW_MTBS_01”.

• layer_id – Layer ID within the service.

• geometry – Geometry for spatial filter. Can be: - A GeoJSON geometry dict (Point, Polygon, Envelope-style bbox) - An Esri JSON geometry string/dict - A bbox string “xmin,ymin,xmax,ymax” - None for no spatial filter

• geometry_type – Esri geometry type. Common values: esriGeometryPoint, esriGeometryEnvelope, esriGeometryPolygon

• spatial_rel – Spatial relationship. Default: esriSpatialRelIntersects.

• where – SQL WHERE clause. Default: “1=1” (all features).

• out_fields – Comma-separated field names or “*” for all. Include the layer’s object-id field (usually “objectid”) if you need stable feature ids: EDW omits the GeoJSON id when that field isn’t requested, leaving only a positional index that shifts between calls.

• max_features – Maximum features to return (capped at server max).

• out_sr – Output spatial reference WKID. Default: 4326 (WGS84).

• return_count_only – If True, return only the count of matching features.

• in_sr – Spatial reference WKID the input geometry is expressed in. Defaults to 4326 (WGS84), matching GeoJSON’s own coordinate convention. This is deliberately independent of out_sr - asking for output in another projection must not change how the input geometry is read.

Returns:

N} if return_count_only.

Return type:

GeoJSON FeatureCollection dict, or {“count”

Note

Some EDW layers 500 or silently return zero features (exceededTransferLimit=True with an empty features list) for a spatial query against a geometrically complex AOI - many disjoint polygon parts and/or holes, e.g. a national forest boundary - even though the AOI itself is valid (returnCountOnly on the identical geometry reports the correct match count). This is an upstream ArcGIS Server limitation, not a client-side bug. When detected, this function automatically falls back to fetching the matching object IDs (which isn’t subject to the same limitation) and then re-fetching features by ID with a plain attribute filter.

Example:

from georest.esri import edw

# One fire's perimeter, as GeoJSON
cameron = edw.query_features("EDW_MTBS_01", 63,
                             where="fire_name = 'CAMERON PEAK'",
                             out_fields="fire_name,year,acres")

# Every MTBS fire inside a bounding box (west, south, east, north)
box = edw.query_features("EDW_MTBS_01", 63, geometry="-106.0,40.4,-105.4,40.9",
                         out_fields="fire_name,year,acres")
len(box["features"])

# Just the count
edw.query_features("EDW_MTBS_01", 63, where="year = 2020", return_count_only=True)

georest.esri.edw.query_features_analytic(service_name: str, layer_id: int, out_analytics: list[dict[str, Any]], partition_by: str | list[str] | None = None, order_by_fields: str | None = None, where: str = '1=1', analytic_where: str | None = None, geometry: dict | str | None = None, geometry_type: str = 'esriGeometryEnvelope', spatial_rel: str = 'esriSpatialRelIntersects', out_fields: str = '*', out_sr: int = 4326, return_geometry: bool = True, in_sr: int | None = None) → dict

Run SQL window-function analytics (rank, running total, LAG/LEAD, …) against an EDW layer via the ArcGIS queryAnalytic operation.

Unlike query_features, rows are not collapsed - each analytic value is appended as a new field on its source feature. Use analytic_where to filter down to specific computed rows (e.g. “rank_val = 1” to keep only the top-ranked feature per partition).

Parameters:

• service_name – Short service name, e.g. “EDW_TimberHarvest_01”.

• layer_id – Layer ID within the service.

• out_analytics –

  List of analytic definitions, each a dict with: - type (str): Esri analytic type, e.g. “RANK”, “SUM”, “LAG”,

  ”ROW_NUMBER”, “PERCENTILE_CONT”.

  • field (str, optional): Input field the analytic operates on (omit for COUNT/ROW_NUMBER).

  • out_name (str): Output field name for the computed value.

  • order_by (str, optional): Field ordering within each partition, e.g. “ACRES DESC” (required by ranking and LAG/LEAD functions).

  • params (dict, optional): Extra analyticParameters, e.g. {“offset”: 1} for LAG/LEAD, {“buckets”: 4} for NTILE.

• partition_by – Field name(s) to group rows into separate windows, e.g. “admin_district_name” or [“admin_forest_name”, “fy_completed”]. Omit to treat the whole result set as one partition.

• order_by_fields – Sort order applied before analytics are computed (separate from an analytic’s own order_by).

• where – SQL WHERE clause filtering the source rows.

• analytic_where – SQL filter applied after analytics are computed, evaluated against output field names (e.g. “acres_rank = 1”).

• geometry – Spatial AOI filter, same as query_features.

• geometry_type – Spatial AOI filter, same as query_features.

• spatial_rel – Spatial AOI filter, same as query_features.

• out_fields – Comma-separated field names or “*” for all.

• out_sr – Output spatial reference WKID. Default: 4326 (WGS84).

• return_geometry – If False, omit geometry from results (cheaper when only the computed analytic values are needed).

• in_sr – Spatial reference WKID of the input geometry, independent of out_sr. Default: 4326 (WGS84).

Returns:

GeoJSON FeatureCollection dict with analytic fields appended to each feature’s properties.

Example:

from georest.esri.edw import query_features_analytic

# Largest timber harvest per ranger district on one forest
# (layer 8 = "Timber Harvest (All Years)")
result = query_features_analytic(
    "EDW_TimberHarvest_01", 8,
    out_analytics=[{
        # RANK takes no input field - the field being ranked is
        # named in order_by.
        "type": "RANK",
        "order_by": "gis_acres DESC", "out_name": "acres_rank",
    }],
    partition_by="admin_district_name",
    analytic_where="acres_rank = 1",
    where="admin_forest_name = 'Medicine Bow-Routt National Forest'"
          " AND gis_acres IS NOT NULL",
    out_fields="sale_name,admin_district_name,gis_acres",
)

georest.esri.edw.top_n_per_group(service_name: str, layer_id: int, field: str, group_by: str | list[str], n: int = 1, descending: bool = True, where: str = '1=1', geometry: dict | str | None = None, geometry_type: str = 'esriGeometryEnvelope', spatial_rel: str = 'esriSpatialRelIntersects', out_fields: str = '*', out_sr: int = 4326, in_sr: int | None = None) → dict

Get the top N features by field within each group_by partition.

Convenience wrapper around query_features_analytic using RANK() + analytic_where to collapse results to just the top N rows per group, e.g. the largest timber harvest per ranger district.

Parameters:

• service_name – Short service name, e.g. “EDW_TimberHarvest_01”.

• layer_id – Layer ID within the service.

• field – Field to rank by, e.g. “gis_acres”.

• group_by – Field name(s) defining each group/partition, e.g. “admin_district_name” or [“admin_forest_name”, “fy_completed”].

• n – How many top features to keep per group. Default: 1.

• descending – If True (default), rank largest-first (top = max value). If False, rank smallest-first (top = min value).

• where – SQL WHERE clause filtering the source rows.

• geometry – Spatial AOI filter, same as query_features.

• geometry_type – Spatial AOI filter, same as query_features.

• spatial_rel – Spatial AOI filter, same as query_features.

• out_fields – Comma-separated field names or “*” for all.

• out_sr – Output spatial reference WKID. Default: 4326 (WGS84).

• in_sr – Spatial reference WKID of the input geometry, independent of out_sr. Default: 4326 (WGS84).

Returns:

GeoJSON FeatureCollection dict containing only the top N features per group, with a rank_val field appended to each.

Note: RANK gives tied values the same rank, so a group with ties at the cutoff returns more than n rows (three features tied at rank 1 all satisfy rank_val <= 1).

Example:

from georest.esri.edw import top_n_per_group

# Largest timber harvest per ranger district on one forest
# (layer 8 = "Timber Harvest (All Years)")
result = top_n_per_group(
    "EDW_TimberHarvest_01", 8,
    field="gis_acres", group_by="admin_district_name",
    where="admin_forest_name = 'Medicine Bow-Routt National Forest'"
          " AND gis_acres IS NOT NULL",
    out_fields="sale_name,admin_district_name,gis_acres",
)

georest.esri.edw.query_features_with_pagination(service_name: str, layer_id: int, geometry: dict | str | None = None, geometry_type: str = 'esriGeometryEnvelope', spatial_rel: str = 'esriSpatialRelIntersects', where: str = '1=1', out_fields: str = '*', max_features: int = 5000, out_sr: int = 4326, in_sr: int | None = None) → dict

Query features with automatic pagination past the server’s per-request cap.

Same args as query_features, but max_features can exceed the layer’s maxRecordCount - results are fetched a page at a time and combined. Returns a combined GeoJSON FeatureCollection.

Note: pages are taken by resultOffset without an explicit sort. ArcGIS doesn’t guarantee a stable row order across requests absent orderByFields, so a layer that reorders between pages could in principle duplicate or skip rows.

Example:

from georest.esri import edw

# More than the server's 2,000-record page: 4,000 MTBS fire locations
# (layer 62, points), fetched a page at a time
fires = edw.query_features_with_pagination("EDW_MTBS_01", 62,
                                           out_fields="fire_name,ig_date,acres",
                                           max_features=4000)
len(fires["features"])   # 4000