imageserver-export 0.2.0

imageserver-export is a small, typed, async Python library that turns an area of any ArcGIS ImageServer into one clean GeoTIFF. Its first use is the open server of the European Environment Agency, which publishes the Copernicus Land Monitoring Service maps of all of Europe, at 10 m: how much of the ground is sealed by concrete and asphalt, how dense the trees are, where the grassland, the water and the wetlands are. No account, no key, no manual download.

async with ImageServer(eea("HRL_ImperviousnessDensity_2018")) as server:
    await server.export(Grid.from_lonlat((8.95, 45.33, 9.45, 45.62), res=100), Path("milan.tif"))

Three dependencies: httpx, numpy, rasterio.

Contents

• What is an ImageServer, in 30 seconds
• What this library does
• Who is it for
• Quick start
• The open maps of Europe
• How it works
• Why nearest neighbour
• Recipes
• For AI agents: the MCP server
• The traps we already fell into
• API at a glance
• Data licence and credits
• Contributing

What is an ImageServer, in 30 seconds

Many public agencies publish their maps with ArcGIS Server. An ImageServer is the part of it that serves rasters: grids of numbers, one per pixel, like satellite-derived land cover. You give it a box, a size and a coordinate system, and its exportImage operation sends back a picture of that box.

That is enough to view a map, but not to use it. One request returns at most a few thousand pixels a side, so Italy at 10 m takes thousands of them. They have to line up exactly, or the pieces shift against each other. And when something goes wrong, the server still answers "200 OK", with an error message where the image should be.

This library does that part for you.

What this library does

Who is it for

• Researchers and analysts who need a Copernicus Land layer for a region, at the resolution of their own grid, in a script they can run again.
• Risk and climate tools that turn land cover into a factor per cell: sealed soil makes rain run off faster, tree cover feeds fires, wetlands change how floods spread.
• Anyone with an ArcGIS ImageServer to read, from the EEA or from a regional or national agency: the library is not tied to one server.

It was written for Limen, an open early-warning platform for landslides, floods and wildfires in Italy, which uses sealed soil to amplify the rain-driven flood trigger.

Quick start

1. Install

pip install imageserver-export

2. Export an area

import asyncio
from pathlib import Path

from imageserver_export import Grid, ImageServer, eea


async def main() -> None:
    async with ImageServer(eea("HRL_ImperviousnessDensity_2018")) as server:
        # Milan and its belt, at 100 m, in LAEA Europe (EPSG:3035).
        grid = Grid.from_lonlat((8.95, 45.33, 9.45, 45.62), res=100)
        result = await server.export(grid, Path("milan.tif"), valid_range=(0, 100))
        print(result.width, "x", result.height, "pixels,", f"{result.valid_fraction:.0%} valid")


asyncio.run(main())

valid_range=(0, 100) tells the library which values are real data: this layer uses 254 for "unclassifiable" and 255 for "no data" (the sea, or outside Europe), and both are kept in the file as they are.

3. Read it

The result is an ordinary GeoTIFF: open it in QGIS, or read it with rasterio:

import rasterio

with rasterio.open("milan.tif") as src:
    band = src.read(1)
    sealed = band[band <= 100]
    print(f"mean sealed soil: {sealed.mean():.1f} %")  # ~25 % around Milan

The open maps of Europe

The EEA serves 26 Copernicus Land layers this way (eea(name) builds the URL, list_services(f"{EEA_URL}/{EEA_FOLDER}") lists them). Among them:

Each layer has a raster attribute table: await server.attribute_table() returns, for the imperviousness layer, 102 rows such as {"Value": 1, "Class_name": "1% imperviousness value", "Area_km2": 222.19}.

How it works

flowchart LR
    A["Grid.from_lonlat(box, res)"] --> B["plan(grid, tile_px)"]
    B --> C{{"exportImage<br/>per tile, 3 at a time"}}
    C -->|"GeoTIFF bytes"| D["check_image()"]
    C -->|"200 + JSON error"| X["ImageServerError"]
    D --> E["Mosaic.paste()<br/>size checked"]
    E --> F["out.tif.part"]
    F -->|"all tiles in"| G["out.tif"]
    F -->|"any failure"| H["deleted"]

1. The grid. Grid.from_lonlat reprojects the box with densified edges (in an equal-area projection the sides of a geographic box curve, and its southern point is not a corner), then grows it outward to whole pixels.
2. The tiles. plan() covers the grid exactly once, edge tiles clamped to it. Every tile's box is a multiple of the resolution, so the server returns exactly the size asked for, with no resampling.
3. The requests. export_params() builds each exportImage query in the grid's own CRS, nearest neighbour by default. check_image() turns a JSON answer into an error with the server's message.
4. The mosaic. Mosaic writes a tiled, deflate-compressed GeoTIFF to .part, refuses a tile whose size differs from the one asked for, counts the valid pixels, writes the provenance tags, and renames the file only on success.

Every step is a plain function or class you can use alone; ImageServer just puts them together.

Why nearest neighbour

Many land-cover layers keep their "no data" marker in the same band as the data: 0 to 100 are percentages, 255 is "nothing here". Bilinear or cubic resampling averages neighbours, so between a 0 % field and a 255 sea pixel it can produce 64: a believable "64 % sealed", exactly where nothing was measured. Nothing downstream can tell it from a real value.

So the default is RSP_NearestNeighbor: every output pixel is one the server really has. Pass interpolation="RSP_BilinearInterpolation" only for continuous layers with no sentinel values.

Recipes

Tree cover of a region, at 20 m

grid = Grid.from_lonlat((15.3, 39.9, 16.9, 41.2), res=20)  # Basilicata
async with ImageServer(eea("HRL_TreeCoverDensity_2018")) as server:
    await server.export(
        grid,
        Path("basilicata_trees.tif"),
        valid_range=(0, 100),
        on_tile=lambda done, total: print(f"{done}/{total}"),
    )

A whole country, without asking for the sea

grid = Grid.from_lonlat((6.6, 35.5, 18.6, 47.1), res=100)  # Italy
async with ImageServer(eea("HRL_ImperviousnessDensity_2018")) as server:
    result = await server.export(
        grid, Path("italy.tif"), valid_range=(0, 100), skip_empty=True, resume=True
    )
print(result.skipped, "of", result.tiles, "tiles were only sea")

skip_empty asks first for one overview, 32 times coarser at the default tile size, and sends no request for a tile with no data on it or next to it: those are filled with "no data". Around Lampedusa it skipped 6 tiles out of 9, with the same pixels as the full export. A feature smaller than one overview pixel and surrounded by sea can be missed, so it is off by default.

resume keeps every tile that arrived in italy.tif.tiles/: if the export fails half way, running it again only asks for the missing ones.

From 100 m pixels to 1 km cells

from imageserver_export import aggregate

cells = aggregate(Path("italy.tif"), 10, valid_range=(0, 100))
cells.mean  # mean sealed soil per 1 km cell, NaN where no data
cells.valid_fraction  # share of each cell that held real data
cells.grid  # the 1 km grid, same corner as the raster

Sentinels stay out of the mean: a coastal cell half in the sea is the mean of its land, with valid_fraction 0.5, not 128 % sealed. The raster is read one strip of cells at a time.

Your own grid, in your CRS

# Snap to an existing 1 km grid in UTM 32N.
grid = Grid.snap((500_000, 4_500_000, 560_000, 4_560_000), res=1000, crs=32632)

Label the classes

async with ImageServer(eea("HRL_WaterWetness_2018")) as server:
    rows = await server.attribute_table()


# The label field is `Class_name` in some layers and `class_name` in others.
def label(row: dict) -> str | None:
    return next((v for k, v in row.items() if k.lower() == "class_name"), None)


labels = {row["Value"]: label(row) for row in rows}
# {0: 'Dry', 1: 'Permanent water', 2: 'Temporary water', 3: 'Permanent wet', ...}

Any other ImageServer

server = ImageServer("https://example.org/arcgis/rest/services/Elevation/ImageServer")
info = await server.info()  # pixel type, size, extent, largest image
await server.export(
    Grid.snap(box, res=info.pixel_size[0], crs=3857),
    Path("dem.tif"),
    pixel_type="F32",
    dtype="float32",
    nodata=None,
)

Bring your own HTTP

import httpx

async with httpx.AsyncClient(timeout=120) as client:  # your timeouts, proxies
    server = ImageServer(url, client, retries=5)


async def send(method: str, url: str, **kw) -> httpx.Response:
    return await my_retrying_client.request(method, url, **kw)  # your policy


server = ImageServer(url, send=send)  # every request goes through `send`

For AI agents: the MCP server

pip install "imageserver-export[mcp]"
claude mcp add imageserver -- imageserver-mcp     # Claude Code

Or in Claude Desktop's claude_desktop_config.json:

{ "mcpServers": { "imageserver": { "command": "imageserver-mcp" } } }

Then ask "how much of the ground is sealed around Milan, and what is the value at the Duomo?". The agent gets six tools:

A service is an EEA layer name or the URL of any ImageServer. Every number comes from the library, not from the model: the agent finds the data and explains it, the same rule as in Limen. Areas are capped (4 million pixels for a summary, 25 million for an export) and a refusal says which resolution would fit, so a single call cannot ask a server for a continent at 10 m.

Over HTTP: imageserver-mcp --transport streamable-http --port 8768 serves /mcp, and /health for container healthchecks. Every option also reads an environment variable (IMAGESERVER_MCP_TRANSPORT, _HOST, _PORT, _ALLOWED_HOSTS, _OUTPUT).

Tried live: the Duomo of Milan is 100% imperviousness value; the box of Milan and its belt at 100 m is 27.9 % sealed on average, 59 % of it not sealed at all.

The traps we already fell into

• "200 OK" is not success. An ImageServer answers a bad request with HTTP 200 and a JSON error. Read as a TIFF, it is just "a corrupt file".
• Four corners are not a box. Reprojecting only the corners of Italy's longitude/latitude box to LAEA cut 6.8 km off the south: Lampedusa ended outside the raster, and 28 grid cells had no data.
• Rounding a box cuts it. Edges are grown to the next pixel, never rounded, or a strip of the area asked for goes missing.
• A tile of the wrong size shifts everything. It is refused, not pasted.
• A half-written file looks whole. Hence .part and the rename.
• Field names are not consistent. The label column of the attribute table is Class_name in the imperviousness layer and class_name in the water and wetness one: compare names case-insensitively.

API at a glance

Full reference: agentengineering.it/imageserver-export.

Data licence and credits

The library is Apache-2.0. The data are not part of it: each server sets its own terms. The Copernicus Land Monitoring Service layers the EEA serves are free and open; credit them as © European Union, Copernicus Land Monitoring Service, European Environment Agency (EEA) when you publish anything made from them.

Contributing

uv sync
make check                                   # ruff, mypy --strict, pytest (100 % coverage)
uv run --group docs python docs/make_images.py  # redraw the README images from live data
uv run --group docs mkdocs serve             # the documentation site, locally

Ideas: a command-line imageserver-export entry point; a bbox given in any CRS; multi-band exports.

License

Apache-2.0.