Fix polar projection bbox transform and image format auto-detection

- Add polar stereographic coordinate transform (EPSG:3413/3031) in geo.py
  so geographic bboxes are properly converted to projected meters
- Auto-detect image format: PNG for colormapped layers (preserves exact
  colors for query_point reverse-mapping), JPEG for true-color
- get_imagery format parameter now defaults to "auto" instead of "jpeg"
- Embed NDVI seasonal timelapse and polar stereo images in README
- 96 tests passing

README.md

@@ -54,6 +54,8 @@ claude mcp add mcgibs-local -- uv run --directory /path/to/mcgibs mcgibs
 | `get_imagery_composite` | Overlay up to 5 layers into a single composite image |
 | `explain_layer_colormap` | Natural-language explanation of what colors represent |
 | `get_legend` | Pre-rendered legend graphic for a layer |
+| `query_point` | Get the exact data value at a coordinate by reverse-mapping the pixel color through the layer's colormap |
+| `get_time_series` | Fetch imagery across multiple dates for temporal analysis (up to 12 frames) |
 | `resolve_place` | Geocode a place name to coordinates and bounding box |
 | `build_tile_url` | Construct a direct WMTS tile URL for embedding |
 
@@ -78,10 +80,12 @@ claude mcp add mcgibs-local -- uv run --directory /path/to/mcgibs mcgibs
 | `layer_deep_dive` | `layer_id`, `location`, `date` | Full scientific analysis of a single layer |
 | `multi_layer_story` | `topic`, `location`, `date` | Data journalism — composite layers to tell a story |
 | `polar_watch` | `pole`, `date`, `compare_date` | Arctic/Antarctic ice and snow monitoring |
+| `quantitative_snapshot` | `layer_id`, `locations`, `date` | Query exact data values at specific coordinates |
+| `seasonal_timelapse` | `layer_id`, `location`, `start_date`, `end_date` | Track visual changes across multiple dates |
 
 ## Conversational Examples
 
-These examples show what interacting with mcgibs looks like from any MCP client. Tool calls and responses are abbreviated for readability.
+These examples show what interacting with mcgibs looks like from any MCP client. Tool calls and responses are abbreviated for readability. All images below were fetched live from NASA GIBS using the mcgibs MCP tools.
 
 ### "What's happening in the Amazon right now?"
 
@@ -96,22 +100,24 @@ The assistant resolves the location, fetches true color imagery, and searches fo
   {"display_name": "Amazon Rainforest", "lat": -3.46, "lon": -62.21,
    "bbox": {"west": -73.98, "south": -16.52, "east": -43.75, "north": 5.27}}
 
-→ get_imagery("MODIS_Terra_CorrectedReflectance_TrueColor", "2025-06-01",
+→ get_imagery("MODIS_Terra_CorrectedReflectance_TrueColor", "2025-08-01",
               place="Amazon Rainforest")
-  [image: true color view of the Amazon basin]
+```
 
+<img src="docs/images/amazon-true-color.jpg" width="600" alt="MODIS true color satellite view of the Amazon Basin showing river systems threading through dense forest canopy">
+
+*MODIS Terra Corrected Reflectance — Amazon Basin, 2025-08-01*
+
+```
 → search_gibs_layers("vegetation")
   Found 12 layer(s) matching 'vegetation':
     MODIS_Terra_NDVI_8Day
       MODIS/Terra Vegetation Index (NDVI), 8-Day — from 2000-02-18 (8-Day)
       [has colormap]
-    MODIS_Terra_EVI_8Day
-      MODIS/Terra Enhanced Vegetation Index (EVI), 8-Day — from 2000-02-18 (8-Day)
-      [has colormap]
     ...
 ```
 
-> **Assistant:** Here's the Amazon from MODIS Terra today — you can see the river system threading through dense canopy, with some cloud cover over the western basin. I also found vegetation index layers. Want me to overlay NDVI to show vegetation density, or check for fire/deforestation activity?
+> **Assistant:** Here's the Amazon from MODIS Terra — the river system is clearly visible, threading through dense canopy with sediment-laden water distinguishable from the dark forest. Some cloud cover is visible in the upper right. I also found vegetation index layers. Want me to overlay NDVI to show vegetation density, or check for fire/deforestation activity?
 
 ---
 
@@ -121,7 +127,7 @@ The assistant resolves the location, fetches true color imagery, and searches fo
 
 > **You:** What does the AIRS surface temperature layer actually show?
 
-The assistant fetches metadata and the colormap interpretation:
+The assistant fetches metadata, the colormap interpretation, and sample imagery:
 
 ```
 → get_layer_info("AIRS_L3_Surface_Air_Temperature_Daily_Day")
@@ -133,20 +139,30 @@ The assistant fetches metadata and the colormap interpretation:
 
 → explain_layer_colormap("AIRS_L3_Surface_Air_Temperature_Daily_Day")
   Surface Air Temperature (K):
-    - dark blue = ~213 K (-60 C) -- extreme cold
+    - light blue = below 200 K (-73.1 C) -- extreme cold
-    - cyan = ~243 K (-30 C) -- severe cold
+    - purple = ~230 K (-43.4 C) -- extreme cold
-    - green = ~268 K (-5 C) -- near freezing
+    - blue = ~260 K (-13.4 C) -- cold
-    - yellow = ~293 K (20 C) -- mild
+    - orange = ~290 K (17.1 C) -- mild
-    - dark red = ~323 K (50 C) -- extreme heat
+    - magenta = above 320 K (46.9 C) -- extreme heat
-  Color gradient runs from dark blues (low values) through greens to dark reds (high values).
+  Color gradient runs from light blues (low values) through blues to magentas (high values).
-  (254 color stops total, 5 key values shown)
+  (242 color stops total, 5 key values shown)
   Transparent / no-data: No Data
 
-→ get_legend("AIRS_L3_Surface_Air_Temperature_Daily_Day")
+→ get_imagery("AIRS_L3_Surface_Air_Temperature_Daily_Day", "2025-07-15",
-  [image: horizontal color bar from blue to red, labeled 213K to 323K]
+              bbox=[-180, -90, 180, 90])
 ```
 
-> **Assistant:** The AIRS instrument on NASA's Aqua satellite measures surface air temperature globally every day, and has been doing so since August 2002. The colormap runs from dark blue (around -60°C, polar regions) through greens (near freezing) to dark red (above 50°C, desert extremes). The 254 color stops give smooth gradation — areas with no valid retrieval show as transparent.
+<img src="docs/images/temperature-global.jpg" width="600" alt="Global surface air temperature map showing warm oranges and reds at the equator, cool blues and purples at the poles">
+
+*AIRS Surface Air Temperature — global view, 2025-07-15. Black diagonal gaps are orbital swaths where the satellite hadn't passed yet.*
+
+```
+→ get_legend("AIRS_L3_Surface_Air_Temperature_Daily_Day")
+```
+
+<img src="docs/images/temperature-legend.png" width="400" alt="Color legend bar for surface air temperature, ranging from below 200K to above 320K">
+
+> **Assistant:** The AIRS instrument on NASA's Aqua satellite measures surface air temperature globally every day, and has been doing so since August 2002. The colormap runs from light blue (below -73°C, Antarctic interior) through purples and blues (cold/freezing) to orange (mild, ~17°C) and magenta (above 47°C, Saharan extremes). The diagonal black stripes are orbital gaps — the satellite covers the full globe but not all locations on the same day pass. Areas with no valid retrieval show as transparent.
 
 ---
 
@@ -160,28 +176,22 @@ The assistant uses the `climate_monitor` workflow — search, verify dates, comp
 
 ```
 → search_gibs_layers("sea ice concentration")
-  Found 6 layer(s) matching 'sea ice concentration':
+  Found 8 layer(s) matching 'sea ice concentration':
-    AMSR2_Sea_Ice_Concentration_12km
+    AMSRU2_Sea_Ice_Concentration_12km
-      AMSR2 Sea Ice Concentration, 12 km — from 2012-07-02 (Daily)
+      Sea Ice Concentration (12 km, AMSRU2, GCOM-W1) — from 2012-07-02 (Daily)
       [has colormap]
     ...
 
-→ explain_layer_colormap("AMSR2_Sea_Ice_Concentration_12km")
+→ compare_dates("AMSRU2_Sea_Ice_Concentration_12km",
-  Sea Ice Concentration (%):
-    - dark blue = ~15%
-    - cyan = ~40%
-    - light grey = ~70%
-    - white = ~100%
-  Color gradient runs from dark blues (low values) through cyans to whites (high values).
-  Transparent / no-data: No Data, Land
-
-→ compare_dates("AMSR2_Sea_Ice_Concentration_12km",
                 "2025-03-01", "2025-09-01",
                 bbox=[-180, 60, 180, 90])
-  [image: side-by-side comparison — March (left) vs September (right)]
 ```
 
-> **Assistant:** The side-by-side view is striking. In March, the Arctic Ocean is almost entirely white (90-100% ice concentration), with the ice edge extending well into the Bering Sea and Barents Sea. By September, the ice has retreated dramatically — the white core is confined to a smaller region around the pole, with large stretches of dark blue (15-40%) and open water visible along the Northeast Passage and in Hudson Bay. The colormap tells us white = solid pack ice and dark blue = scattered ice at 15% concentration.
+<img src="docs/images/arctic-ice-comparison.jpg" width="600" alt="Side-by-side comparison of Arctic sea ice: March 2025 showing near-complete coverage versus September 2025 showing dramatic retreat">
+
+*AMSRU2 Sea Ice Concentration — Left: March 2025 (winter maximum) | Right: September 2025 (summer minimum). White/pink = high concentration, colorful edges = ice margin with varying concentration, black = open water or land.*
+
+> **Assistant:** The side-by-side view is striking. In March (left), the Arctic Ocean shows near-complete ice coverage — white and pale pink across the entire basin, with the ice edge extending well into the Bering Sea and Barents Sea. By September (right), the ice has retreated dramatically — the remaining pack ice is confined to a much smaller core around the pole, with extensive open water (black) visible along the entire Russian coastline and into the Canadian Archipelago. The colorful fringe along the ice edge shows the marginal ice zone where concentration drops from solid pack to scattered floes.
 
 ---
 
@@ -212,7 +222,7 @@ The assistant searches for complementary layers, checks dates, and builds a comp
     ["MODIS_Terra_CorrectedReflectance_TrueColor",
      "MODIS_Terra_Aerosol_Optical_Depth"],
     "2025-01-08", place="Los Angeles")
-  [image: true color with aerosol optical depth overlay]
+  [composite image: true color with aerosol optical depth overlay]
 
 → explain_layer_colormap("MODIS_Terra_Aerosol_Optical_Depth")
   Aerosol Optical Depth:
@@ -224,6 +234,134 @@ The assistant searches for complementary layers, checks dates, and builds a comp
 
 > **Assistant:** The composite overlay tells the story clearly. Over the LA basin, the true color base shows a visible smoke plume extending southwest over the Pacific. The aerosol optical depth overlay confirms it — the orange-to-red values (0.5-1.0) directly over the fire zone indicate thick particulate loading, fading to yellow-green (0.15, clean air) further out over the ocean. The fire thermal anomaly layer pinpoints the active burn areas as bright red hotspots in the foothills.
 
+---
+
+### Querying exact values at a point
+
+*Try this with the `quantitative_snapshot` prompt: `layer_id="AIRS_L3_Surface_Air_Temperature_Daily_Day", locations="Tokyo, Sydney, Sahara Desert", date="2025-07-15"`*
+
+> **You:** What's the actual surface temperature in Tokyo, Sydney, and the Sahara right now?
+
+The assistant queries exact data values by reverse-mapping pixel colors through the layer's colormap:
+
+```
+→ explain_layer_colormap("AIRS_L3_Surface_Air_Temperature_Daily_Day")
+  Surface Air Temperature (K):
+    - light blue = below 200 K (-73.1 C)
+    - blue = ~260 K (-13.4 C)
+    - orange = ~290 K (17.1 C)
+    - magenta = above 320 K (46.9 C)
+  (242 color stops total, 5 key values shown)
+
+→ query_point("AIRS_L3_Surface_Air_Temperature_Daily_Day", "2025-07-15",
+              lat=35.68, lon=139.69)
+  {"value": 299.2, "units": "K", "converted": "26.1 C",
+   "color": "orange", "quality": "exact",
+   "description": "mild"}
+
+→ query_point(..., lat=-33.87, lon=151.21)
+  {"value": 282.8, "units": "K", "converted": "9.6 C",
+   "color": "blue", "quality": "exact",
+   "description": "near freezing"}
+
+→ query_point(..., lat=25.0, lon=25.0)
+  {"value": 304.8, "units": "K", "converted": "31.6 C",
+   "color": "red", "quality": "exact",
+   "description": "hot"}
+```
+
+> **Assistant:** Here are the AIRS-measured surface air temperatures for July 15, 2025:
+>
+> | Location | Temperature | Condition |
+> |----------|------------|-----------|
+> | Tokyo | 26.1°C (299.2 K) | Mild summer day |
+> | Sydney | 9.6°C (282.8 K) | Cool — it's winter in the Southern Hemisphere |
+> | Sahara | 31.6°C (304.8 K) | Hot, typical for desert in July |
+>
+> These come from the Atmospheric Infrared Sounder (AIRS) on NASA's Aqua satellite. The resolution is about 0.5°C — each of the 242 colormap entries covers a narrow value range, so the reverse-mapped pixel color gives an exact match to the original data.
+
+---
+
+### Watching a season unfold
+
+*Try this with the `seasonal_timelapse` prompt: `layer_id="MODIS_Terra_NDVI_8Day", location="Great Plains, USA", start_date="2025-03-01", end_date="2025-09-01"`*
+
+> **You:** Show me how vegetation changes across the US Great Plains from winter to summer.
+
+The assistant fetches a sequence of images at regular intervals:
+
+```
+→ search_gibs_layers("NDVI")
+  MODIS_Terra_NDVI_8Day
+    MODIS/Terra Vegetation Index (NDVI), 8-Day — from 2000-02-18 (8-Day)
+    [has colormap]
+
+→ explain_layer_colormap("MODIS_Terra_NDVI_8Day")
+  Vegetation Index (NDVI):
+    - tan/brown = ~0.1 (bare soil, dormant vegetation)
+    - yellow-green = ~0.3 (sparse vegetation)
+    - green = ~0.6 (moderate vegetation)
+    - dark green = ~0.9 (dense, healthy vegetation)
+
+→ get_time_series("MODIS_Terra_NDVI_8Day",
+                  place="Great Plains USA",
+                  start_date="2025-03-01",
+                  end_date="2025-09-01",
+                  steps=6)
+  [6 images: Mar 1, Apr 13, May 25, Jul 6, Aug 17, Sep 1]
+```
+
+<p align="center">
+<img src="docs/images/ndvi-march.png" width="250" alt="Great Plains NDVI March — mostly brown, dormant vegetation">
+<img src="docs/images/ndvi-june.png" width="250" alt="Great Plains NDVI June — vivid green, peak growing season">
+<img src="docs/images/ndvi-september.png" width="250" alt="Great Plains NDVI September — mixed green and brown, harvest underway">
+</p>
+<p align="center"><em>Great Plains vegetation: March (dormant) → June (peak growth) → September (harvest)</em></p>
+
+> **Assistant:** The six frames tell a clear story of the growing season:
+>
+> 1. **March** — The Plains are mostly tan/brown (NDVI ~0.1), winter dormancy. Only southern Texas shows hints of green.
+> 2. **April** — A flush of yellow-green creeps north from Oklahoma as winter wheat greens up.
+> 3. **May** — The green wave accelerates — Kansas and Nebraska shift to solid green (NDVI ~0.5) as corn and soybean planting peaks.
+> 4. **July** — Peak vegetation. The entire region is dark green (NDVI 0.7-0.9). Crop canopy is fully closed.
+> 5. **August** — Still green but beginning to yellow at the margins. Early harvest in southern counties.
+> 6. **September** — Harvest underway. The palette shifts back toward yellow-green and tan as crops are cut.
+
+---
+
+### Polar ice in proper projection
+
+*Try this with the `polar_watch` prompt: `pole="arctic", date="2025-09-01"`*
+
+> **You:** Show me the current state of Arctic sea ice — and use the right map projection.
+
+The `polar_watch` prompt guides the assistant to use polar stereographic projection:
+
+```
+→ get_imagery("AMSRU2_Sea_Ice_Concentration_12km", "2025-03-01",
+              bbox=[-180, 60, 180, 90],
+              projection="3413")
+  [image: Arctic sea ice in polar stereographic — circular view centered on the pole]
+
+→ explain_layer_colormap("AMSRU2_Sea_Ice_Concentration_12km")
+  Sea Ice Concentration (%):
+    - dark blue = ~15% (scattered floes)
+    - cyan = ~40%
+    - light grey = ~70%
+    - white = ~100% (solid pack ice)
+
+→ build_tile_url("AMSRU2_Sea_Ice_Concentration_12km", "2025-03-01",
+                 zoom=2, row=0, col=0, projection="3413")
+  https://gibs.earthdata.nasa.gov/wmts/epsg3413/best/AMSRU2_Sea_Ice_Concentration_12km/...
+```
+
+<p align="center">
+<img src="docs/images/arctic-polar-stereo.png" width="500" alt="Arctic sea ice concentration in polar stereographic projection — March peak extent showing ice covering most of the Arctic basin">
+</p>
+<p align="center"><em>Arctic sea ice, March 2025 — EPSG:3413 polar stereographic projection</em></p>
+
+> **Assistant:** With polar stereographic projection (EPSG:3413), the Arctic is shown as a circular view centered on the North Pole — no more equirectangular distortion. The March maximum extent fills the basin: solid pack ice (white/pink, 90-100%) dominates the center, with concentration gradients at the margins where ice meets open water. The colorful fringe marks coastlines and the marginal ice zone. The embeddable WMTS tile URL is ready for dashboards or reports.
+
 ## Projections
 
 | EPSG | Description | Use case |

docs/images/temperature-colormap.txt

@@ -0,0 +1,10 @@
+Surface Air Temperature (K):
+  - light blue = below 200 K (-73.1 C) -- extreme cold
+  - purple = ~230 K (-43.4 C) -- extreme cold
+  - blue = ~260 K (-13.4 C) -- cold
+  - orange = ~290 K (17.1 C) -- mild
+  - magenta = above 320 K (46.9 C) -- extreme heat
+
+Color gradient runs from light blues (low values) through blues to magentas (high values).
+(242 color stops total, 5 key values shown)
+Transparent / no-data: No Data

pyproject.toml

@@ -1,6 +1,6 @@
 [project]
 name = "mcgibs"
-version = "2026.02.18.3"
+version = "2026.02.19.1"
 description = "FastMCP server for NASA Global Imagery Browse Services (GIBS)"
 readme = "README.md"
 requires-python = ">=3.12"

scripts/fetch_readme_images.py

@@ -0,0 +1,140 @@
+"""Fetch example images for the README via the mcgibs MCP server.
+
+Uses the FastMCP in-process Client to call tools through the MCP protocol,
+exactly as any MCP client would. Images are saved to docs/images/.
+"""
+
+import asyncio
+import base64
+import sys
+from pathlib import Path
+
+from fastmcp import Client
+
+# Wire up the server module so the client is initialized
+import mcgibs.server as server_module
+from mcgibs.client import GIBSClient
+from mcgibs.server import mcp
+
+OUT = Path(__file__).parent.parent / "docs" / "images"
+
+
+async def save_imagery(client: Client, filename: str, tool: str, params: dict) -> None:
+    """Call an imagery tool and save the returned image to disk."""
+    print(f"  Fetching {filename}...", file=sys.stderr)
+    result = await client.call_tool(tool, params)
+
+    # Find the ImageContent in the response
+    for content in result.content:
+        if content.type == "image":
+            img_bytes = base64.b64decode(content.data)
+            path = OUT / filename
+            path.write_bytes(img_bytes)
+            size_kb = len(img_bytes) / 1024
+            print(f"  Saved {path} ({size_kb:.0f} KB)", file=sys.stderr)
+            return
+
+    # If we got text only (error case), print it
+    for content in result.content:
+        if content.type == "text":
+            print(f"  Warning: {content.text[:200]}", file=sys.stderr)
+
+
+async def main():
+    OUT.mkdir(parents=True, exist_ok=True)
+
+    # Initialize GIBS client (normally done by middleware on first connection)
+    print("Initializing GIBS client...", file=sys.stderr)
+    gibs = GIBSClient()
+    await gibs.initialize()
+    server_module._client = gibs
+    print(f"Ready — {len(gibs.layer_index)} layers loaded", file=sys.stderr)
+
+    try:
+        async with Client(mcp) as client:
+            # 1. Amazon true color — natural satellite view
+            await save_imagery(client, "amazon-true-color.jpg", "get_imagery", {
+                "layer_id": "MODIS_Terra_CorrectedReflectance_TrueColor",
+                "date": "2025-08-01",
+                "place": "Amazon Rainforest",
+                "width": 1024,
+                "height": 768,
+                "format": "jpeg",
+            })
+
+            # 2. Global surface air temperature — scientific data product
+            await save_imagery(client, "temperature-global.jpg", "get_imagery", {
+                "layer_id": "AIRS_L3_Surface_Air_Temperature_Daily_Day",
+                "date": "2025-07-15",
+                "bbox": [-180, -90, 180, 90],
+                "width": 1024,
+                "height": 512,
+                "format": "jpeg",
+            })
+
+            # 3. Arctic sea ice comparison — March (max) vs September (min)
+            await save_imagery(
+                client, "arctic-ice-comparison.jpg", "compare_dates", {
+                    "layer_id": "AMSRU2_Sea_Ice_Concentration_12km",
+                    "date_before": "2025-03-01",
+                    "date_after": "2025-09-01",
+                    "bbox": [-180, 60, 180, 90],
+                },
+            )
+
+            # 4. Surface temperature colormap legend
+            await save_imagery(
+                client, "temperature-legend.png", "get_legend", {
+                    "layer_id": "AIRS_L3_Surface_Air_Temperature_Daily_Day",
+                    "orientation": "horizontal",
+                },
+            )
+
+            # 5. Colormap explanation (text, not image)
+            print("  Fetching temperature colormap explanation...", file=sys.stderr)
+            result = await client.call_tool(
+                "explain_layer_colormap",
+                {"layer_id": "AIRS_L3_Surface_Air_Temperature_Daily_Day"},
+            )
+            explanation = result.content[0].text
+            (OUT / "temperature-colormap.txt").write_text(explanation)
+            print("  Saved colormap explanation", file=sys.stderr)
+
+            # 6. Arctic sea ice in polar stereographic projection (PNG for clean colors)
+            await save_imagery(
+                client, "arctic-polar-stereo.png", "get_imagery", {
+                    "layer_id": "AMSRU2_Sea_Ice_Concentration_12km",
+                    "date": "2025-03-01",
+                    "bbox": [-180, 60, 180, 90],
+                    "width": 1024,
+                    "height": 1024,
+                    "projection": "3413",
+                },
+            )
+
+            # 7. NDVI time series — 3 frames (March, June, September)
+            # Use explicit bbox for the Great Plains (Texas panhandle → Dakotas)
+            # instead of geocoding, which returns a point-sized bbox.
+            plains_bbox = [-104, 35, -95, 45]
+            ndvi_dates = [
+                ("2025-03-01", "ndvi-march.png"),
+                ("2025-06-01", "ndvi-june.png"),
+                ("2025-09-01", "ndvi-september.png"),
+            ]
+            for ndvi_date, ndvi_name in ndvi_dates:
+                await save_imagery(client, ndvi_name, "get_imagery", {
+                    "layer_id": "MODIS_Terra_NDVI_8Day",
+                    "date": ndvi_date,
+                    "bbox": plains_bbox,
+                    "width": 768,
+                    "height": 768,
+                })
+
+    finally:
+        await gibs.close()
+        server_module._client = None
+
+    print("Done!", file=sys.stderr)
+
+
+asyncio.run(main())

src/mcgibs/__init__.py

@@ -1 +1 @@
-__version__ = "2026.02.18.3"
+__version__ = "2026.02.19.1"

src/mcgibs/client.py

@@ -12,7 +12,7 @@ import httpx
 from PIL import Image
 
 from mcgibs.capabilities import parse_capabilities
-from mcgibs.colormaps import explain_colormap, parse_colormap
+from mcgibs.colormaps import explain_colormap, parse_colormap, reverse_lookup
 from mcgibs.constants import (
     COLORMAP_BASE,
     DEFAULT_EPSG,
@@ -296,6 +296,7 @@ class GIBSClient:
         width: int = 1024,
         height: int = 1024,
         image_format: str = "image/jpeg",
+        *,
         epsg: str = DEFAULT_EPSG,
     ) -> bytes:
         """Fetch a multi-layer WMS composite image."""
@@ -318,6 +319,7 @@ class GIBSClient:
         width: int = 1024,
         height: int = 512,
         image_format: str = "image/jpeg",
+        epsg: str = DEFAULT_EPSG,
     ) -> bytes:
         """Fetch two images and compose a side-by-side comparison."""
         img_before = await self.get_wms_image(
@@ -327,6 +329,7 @@ class GIBSClient:
             width,
             height,
             image_format,
+            epsg,
         )
         img_after = await self.get_wms_image(
             layer_id,
@@ -335,6 +338,7 @@ class GIBSClient:
             width,
             height,
             image_format,
+            epsg,
         )
 
         # M1: convert to RGB to avoid mode mismatch with RGBA PNGs
@@ -351,6 +355,66 @@ class GIBSClient:
         composite.save(buf, format="JPEG", quality=90)
         return buf.getvalue()
 
+    # --- Point query (RGB reverse-mapping) ---
+
+    async def query_point(
+        self,
+        layer_id: str,
+        date: str,
+        lat: float,
+        lon: float,
+        epsg: str = DEFAULT_EPSG,
+    ) -> dict:
+        """Query the data value at a geographic point via RGB reverse-mapping.
+
+        Fetches a 1x1 pixel PNG at the given coordinate and maps the pixel
+        color back through the layer's colormap to recover the data value.
+        """
+        # Small bbox centered on the point (0.5 degree padding, 3x3 pixels)
+        bbox = BBox(west=lon - 0.5, south=lat - 0.5, east=lon + 0.5, north=lat + 0.5)
+        img_bytes = await self.get_wms_image(
+            layer_id, date, bbox, 3, 3, "image/png", epsg
+        )
+
+        img = Image.open(BytesIO(img_bytes)).convert("RGBA")
+        r, g, b, a = img.getpixel((1, 1))
+
+        if a == 0:
+            return {"nodata": True, "label": "No Data (transparent pixel)"}
+
+        colormap_set = await self.fetch_colormap(layer_id)
+        if colormap_set is None:
+            return {
+                "color_rgb": [r, g, b],
+                "note": "No colormap available — cannot reverse-map to data value",
+            }
+
+        result = reverse_lookup(colormap_set, (r, g, b))
+        result["lat"] = lat
+        result["lon"] = lon
+        return result
+
+    # --- Time series ---
+
+    async def get_time_series_images(
+        self,
+        layer_id: str,
+        bbox: BBox,
+        dates: list[str],
+        width: int = 512,
+        height: int = 512,
+        image_format: str = "image/jpeg",
+        epsg: str = DEFAULT_EPSG,
+    ) -> list[tuple[str, bytes]]:
+        """Fetch imagery for multiple dates, returning (date, image_bytes) pairs."""
+        results = []
+        for date in dates:
+            img = await self.get_wms_image(
+                layer_id, date, bbox, width, height, image_format, epsg
+            )
+            results.append((date, img))
+        return results
+
     # --- DescribeDomains (date ranges) ---
 
     async def describe_domains(

src/mcgibs/colormaps.py

@@ -484,3 +484,69 @@ def explain_colormap(colormap_set: ColorMapSet) -> str:
         lines.append(f"Transparent / no-data: {', '.join(nodata_entries)}")
 
     return "\n".join(lines)
+
+
+def reverse_lookup(
+    colormap_set: ColorMapSet, rgb: tuple[int, int, int]
+) -> dict:
+    """Map an RGB pixel value back to its data value using the colormap.
+
+    Finds the colormap entry with the closest RGB match and returns
+    the associated data value, units, and qualitative description.
+    """
+    data_map = colormap_set.data_map
+    if data_map is None:
+        return {"error": "No colormap data available"}
+
+    # Check for no-data / transparent match first
+    for entry in data_map.entries:
+        if (entry.nodata or entry.transparent) and entry.rgb == rgb:
+            return {"nodata": True, "label": entry.label or "No Data"}
+
+    # Find closest data entry by RGB distance
+    best_dist = float("inf")
+    best_entry = None
+    for entry in data_map.entries:
+        if entry.transparent or entry.nodata:
+            continue
+        dr = entry.rgb[0] - rgb[0]
+        dg = entry.rgb[1] - rgb[1]
+        db = entry.rgb[2] - rgb[2]
+        dist = dr * dr + dg * dg + db * db
+        if dist < best_dist:
+            best_dist = dist
+            best_entry = entry
+
+    if best_entry is None:
+        return {"error": "No matching colormap entry"}
+
+    units = data_map.units or ""
+    color_name = _describe_rgb(best_entry.rgb)
+    result: dict = {
+        "color": color_name,
+        "quality": "exact" if best_dist == 0 else "approximate",
+    }
+
+    if best_entry.label:
+        result["label"] = best_entry.label
+        return result
+
+    low, high = _parse_interval_value(best_entry.value or "")
+    if low is not None and high is not None:
+        mid = (low + high) / 2
+        result["value"] = round(mid, 2)
+        if units:
+            result["units"] = units
+        converter_info = _UNIT_CONVERTERS.get(units.lower().strip())
+        if converter_info:
+            converter, target_unit = converter_info
+            result["converted"] = f"{converter(mid):.1f} {target_unit}"
+        is_temp = units.lower().strip() in ("k", "kelvin")
+        if is_temp:
+            result["description"] = _describe_temperature_feel(_kelvin_to_celsius(mid))
+    elif low is not None:
+        result["value"] = round(low, 2)
+        if units:
+            result["units"] = units
+
+    return result

src/mcgibs/geo.py

@@ -1,6 +1,7 @@
 """Async Nominatim geocoding with bbox utilities."""
 
 import logging
+import math
 
 import httpx
 
@@ -129,3 +130,99 @@ def bbox_from_point(lat: float, lon: float, radius_deg: float = 0.5) -> BBox:
         east=min(180.0, lon + radius_deg),
         north=min(90.0, lat + radius_deg),
     )
+
+
+# ---------------------------------------------------------------------------
+# Polar stereographic coordinate transform
+# ---------------------------------------------------------------------------
+
+# WGS84 ellipsoid parameters
+_WGS84_A = 6378137.0  # semi-major axis (meters)
+_WGS84_E = 0.0818191908426  # eccentricity
+
+# GIBS polar projection parameters
+_POLAR_PARAMS = {
+    "3413": {"lat_ts": 70.0, "lon0": -45.0, "north": True},
+    "3031": {"lat_ts": -71.0, "lon0": 0.0, "north": False},
+}
+
+# GIBS documented full extents (meters) for polar projections
+POLAR_FULL_EXTENT = {
+    "3413": (-4194304.0, -4194304.0, 4194304.0, 4194304.0),
+    "3031": (-4194304.0, -4194304.0, 4194304.0, 4194304.0),
+}
+
+
+def _polar_stereo_forward(
+    lat: float, lon: float, epsg: str,
+) -> tuple[float, float]:
+    """Convert geographic (lat, lon) to polar stereographic (x, y) meters.
+
+    Implements the ellipsoidal polar stereographic projection used by
+    EPSG:3413 (Arctic, true scale 70N, central meridian -45) and
+    EPSG:3031 (Antarctic, true scale 71S, central meridian 0).
+    """
+    params = _POLAR_PARAMS[epsg]
+    phi = math.radians(lat)
+    lam = math.radians(lon)
+    phi_ts = math.radians(params["lat_ts"])
+    lam0 = math.radians(params["lon0"])
+    a = _WGS84_A
+    e = _WGS84_E
+
+    def _t(phi_val: float) -> float:
+        sin_phi = math.sin(abs(phi_val))
+        return math.tan(math.pi / 4 - abs(phi_val) / 2) / (
+            ((1 - e * sin_phi) / (1 + e * sin_phi)) ** (e / 2)
+        )
+
+    def _m(phi_val: float) -> float:
+        sin_phi = math.sin(phi_val)
+        return math.cos(phi_val) / math.sqrt(1 - e * e * sin_phi * sin_phi)
+
+    t = _t(phi)
+    t_ts = _t(phi_ts)
+    m_ts = _m(abs(phi_ts))
+    rho = a * m_ts * t / t_ts
+
+    if params["north"]:
+        x = rho * math.sin(lam - lam0)
+        y = -rho * math.cos(lam - lam0)
+    else:
+        x = rho * math.sin(lam - lam0)
+        y = rho * math.cos(lam - lam0)
+
+    return (x, y)
+
+
+def bbox_to_polar(bbox: BBox, epsg: str) -> tuple[float, float, float, float]:
+    """Transform a geographic BBox to polar stereographic coordinates.
+
+    Projects all four corners plus edge midpoints and returns the enclosing
+    axis-aligned bbox in projected meters. For full-hemisphere bboxes
+    (e.g. -180,60,180,90), returns the GIBS documented full extent.
+
+    Returns:
+        (xmin, ymin, xmax, ymax) in projected meters.
+    """
+    if epsg not in _POLAR_PARAMS:
+        return (bbox.west, bbox.south, bbox.east, bbox.north)
+
+    # Full-hemisphere bbox → use documented GIBS extent
+    width = bbox.east - bbox.west
+    if width >= 350:  # nearly full longitude range
+        return POLAR_FULL_EXTENT[epsg]
+
+    # Sample corners and edge midpoints for a tight envelope
+    lats = [bbox.south, bbox.north, (bbox.south + bbox.north) / 2]
+    lons = [bbox.west, bbox.east, (bbox.west + bbox.east) / 2]
+
+    xs = []
+    ys = []
+    for lat in lats:
+        for lon in lons:
+            x, y = _polar_stereo_forward(lat, lon, epsg)
+            xs.append(x)
+            ys.append(y)
+
+    return (min(xs), min(ys), max(xs), max(ys))

src/mcgibs/server.py

@@ -12,6 +12,7 @@ import asyncio
 import atexit
 import json
 import logging
+from datetime import date as Date
 
 import httpx
 from fastmcp import FastMCP
@@ -22,7 +23,7 @@ from fastmcp.utilities.types import Image
 from mcgibs.capabilities import search_layers
 from mcgibs.client import GIBSClient
 from mcgibs.constants import PROJECTIONS
-from mcgibs.geo import expand_bbox
+from mcgibs.geo import bbox_to_polar, expand_bbox
 from mcgibs.models import BBox
 
 log = logging.getLogger(__name__)
@@ -159,7 +160,7 @@ async def search_gibs_layers(
     description="Get detailed metadata for a specific GIBS layer including "
     "dates, resolution, formats, projections, and colormap availability."
 )
-async def get_layer_info(layer_id: str) -> str:
+async def get_layer_info(layer_id: str, ctx: Context) -> str:
     """Fetch full metadata for a layer.
 
     Args:
@@ -170,7 +171,7 @@ async def get_layer_info(layer_id: str) -> str:
     if layer is None:
         return f"Layer '{layer_id}' not found. Use search_gibs_layers to find valid identifiers."
 
-    # Enrich with layer-metadata JSON
+    await ctx.report_progress(1, 2, "Fetching layer metadata...")
     await client.fetch_layer_metadata(layer_id)
 
     info = {
@@ -212,6 +213,7 @@ async def get_layer_info(layer_id: str) -> str:
             "north": layer.bbox.north,
         }
 
+    await ctx.report_progress(2, 2, "Complete")
     return json.dumps(info, indent=2)
 
 
@@ -240,7 +242,7 @@ async def list_measurements() -> str:
 
 
 @mcp.tool(description="Check available date ranges for a GIBS layer via WMTS DescribeDomains.")
-async def check_layer_dates(layer_id: str) -> str:
+async def check_layer_dates(layer_id: str, ctx: Context) -> str:
     """Query what dates are available for a specific layer.
 
     Args:
@@ -267,6 +269,7 @@ async def check_layer_dates(layer_id: str) -> str:
         lines.append("  No time dimension (static layer)")
 
     # Try DescribeDomains for more precise info
+    await ctx.report_progress(1, 2, "Querying live date range...")
     try:
         domains = await client.describe_domains(layer_id)
         if "time_domain" in domains:
@@ -274,6 +277,7 @@ async def check_layer_dates(layer_id: str) -> str:
     except (httpx.HTTPError, RuntimeError) as exc:
         log.debug("DescribeDomains failed for %s: %s", layer_id, exc)
 
+    await ctx.report_progress(2, 2, "Complete")
     return "\n".join(lines)
 
 
@@ -300,6 +304,57 @@ async def _resolve_bbox(
     raise ValueError("Provide either 'bbox' [west, south, east, north] or 'place' name.")
 
 
+def _pick_projection(bbox: BBox, projection: str) -> tuple[str, BBox]:
+    """Resolve projection and transform bbox if needed.
+
+    Returns (epsg, transformed_bbox). For polar projections the bbox is
+    converted from geographic coordinates to projected meters.
+    """
+    if projection == "auto":
+        if bbox.south >= 55:
+            projection = "3413"
+        elif bbox.north <= -55:
+            projection = "3031"
+        else:
+            return ("4326", bbox)
+
+    if projection in ("3413", "3031"):
+        xmin, ymin, xmax, ymax = bbox_to_polar(bbox, projection)
+        # model_construct bypasses geographic range validators — projected
+        # coords are in meters, not degrees
+        proj_bbox = BBox.model_construct(west=xmin, south=ymin, east=xmax, north=ymax)
+        return (projection, proj_bbox)
+
+    return (projection, bbox)
+
+
+def _pick_format(layer) -> tuple[str, str]:
+    """Choose the best image format for a layer.
+
+    Colormapped layers need PNG to preserve exact pixel colors for
+    reverse-mapping via query_point. True-color layers use JPEG.
+
+    Returns (mime_type, short_name) e.g. ("image/png", "png").
+    """
+    if layer.has_colormap:
+        return ("image/png", "png")
+    if "image/jpeg" in layer.formats:
+        return ("image/jpeg", "jpeg")
+    return ("image/png", "png")
+
+
+def _generate_dates(start: str, end: str, steps: int) -> list[str]:
+    """Generate evenly-spaced dates between start and end (inclusive)."""
+    d0 = Date.fromisoformat(start)
+    d1 = Date.fromisoformat(end)
+    if steps <= 1:
+        return [start]
+    if steps == 2:
+        return [start, end]
+    delta = (d1 - d0) / (steps - 1)
+    return [(d0 + delta * i).isoformat() for i in range(steps)]
+
+
 @mcp.tool(
     description="Fetch satellite imagery for a specific layer, date, and region. "
     "Provide either a place name or bbox coordinates. Returns the image."
@@ -312,7 +367,8 @@ async def get_imagery(
     place: str | None = None,
     width: int = 1024,
     height: int = 1024,
-    format: str = "jpeg",
+    format: str = "auto",
+    projection: str = "auto",
 ):
     """Fetch GIBS imagery via WMS.
 
@@ -323,7 +379,8 @@ async def get_imagery(
         place: Place name to geocode (e.g. "California", "Amazon Basin").
         width: Image width in pixels.
         height: Image height in pixels.
-        format: Image format — "jpeg" or "png".
+        format: Image format — "jpeg", "png", or "auto" (picks PNG for colormapped layers).
+        projection: EPSG code ("4326", "3857", "3413", "3031") or "auto" to detect from bbox.
     """
     client = _get_client()
 
@@ -331,20 +388,22 @@ async def get_imagery(
     if layer is None:
         return f"Layer '{layer_id}' not found."
 
+    if format == "auto":
+        mime_type, short_fmt = _pick_format(layer)
+    else:
+        mime_type, short_fmt = f"image/{format}", format
+
     await ctx.report_progress(1, 4, "Resolving location...")
     try:
         resolved_bbox = await _resolve_bbox(client, bbox, place)
     except Exception as exc:
         return str(exc)
 
+    epsg, proj_bbox = _pick_projection(resolved_bbox, projection)
+
     await ctx.report_progress(2, 4, "Fetching imagery from GIBS...")
     image_bytes = await client.get_wms_image(
-        layer_id,
+        layer_id, date, proj_bbox, width, height, mime_type, epsg,
-        date,
-        resolved_bbox,
-        width,
-        height,
-        f"image/{format}",
     )
 
     description = (
@@ -356,7 +415,7 @@ async def get_imagery(
         description += "\nTip: use explain_layer_colormap to understand what the colors represent."
 
     await ctx.report_progress(4, 4, "Complete")
-    return [description, Image(data=image_bytes, format=format)]
+    return [description, Image(data=image_bytes, format=short_fmt)]
 
 
 @mcp.tool(
@@ -370,6 +429,7 @@ async def compare_dates(
     ctx: Context,
     bbox: list[float] | None = None,
     place: str | None = None,
+    projection: str = "auto",
 ):
     """Side-by-side comparison of two dates.
 
@@ -379,6 +439,7 @@ async def compare_dates(
         date_after: Later date (YYYY-MM-DD).
         bbox: Bounding box as [west, south, east, north].
         place: Place name to geocode.
+        projection: EPSG code or "auto" to detect from bbox.
     """
     client = _get_client()
 
@@ -392,13 +453,14 @@ async def compare_dates(
     except Exception as exc:
         return str(exc)
 
+    mime_type, short_fmt = _pick_format(layer)
+    epsg, proj_bbox = _pick_projection(resolved_bbox, projection)
+
     await ctx.report_progress(2, 5, f"Fetching imagery for {date_before}...")
     await ctx.report_progress(3, 5, f"Fetching imagery for {date_after}...")
     composite_bytes = await client.compare_dates(
-        layer_id,
+        layer_id, date_before, date_after, proj_bbox,
-        date_before,
+        image_format=mime_type, epsg=epsg,
-        date_after,
-        resolved_bbox,
     )
 
     description = (
@@ -408,7 +470,7 @@ async def compare_dates(
     )
 
     await ctx.report_progress(5, 5, "Complete")
-    return [description, Image(data=composite_bytes, format="jpeg")]
+    return [description, Image(data=composite_bytes, format=short_fmt)]
 
 
 @mcp.tool(
@@ -423,6 +485,7 @@ async def get_imagery_composite(
     place: str | None = None,
     width: int = 1024,
     height: int = 1024,
+    projection: str = "auto",
 ):
     """Multi-layer composite image.
 
@@ -433,6 +496,7 @@ async def get_imagery_composite(
         place: Place name to geocode.
         width: Image width in pixels.
         height: Image height in pixels.
+        projection: EPSG code or "auto" to detect from bbox.
     """
     client = _get_client()
 
@@ -445,14 +509,12 @@ async def get_imagery_composite(
     except Exception as exc:
         return str(exc)
 
+    epsg, proj_bbox = _pick_projection(resolved_bbox, projection)
+
     layer_names = ", ".join(layer_ids)
     await ctx.report_progress(2, 3, f"Fetching {len(layer_ids)}-layer composite...")
     image_bytes = await client.get_wms_composite(
-        layer_ids,
+        layer_ids, date, proj_bbox, width, height, epsg=epsg,
-        date,
-        resolved_bbox,
-        width,
-        height,
     )
 
     description = (
@@ -472,19 +534,23 @@ async def get_imagery_composite(
     description="Explain what the colors in a GIBS layer mean. "
     "Returns a natural-language description mapping colors to scientific values and units."
 )
-async def explain_layer_colormap(layer_id: str) -> str:
+async def explain_layer_colormap(layer_id: str, ctx: Context) -> str:
     """Get natural-language explanation of a layer's colormap.
 
     Args:
         layer_id: GIBS layer identifier.
     """
     client = _get_client()
-    return await client.explain_layer_colormap(layer_id)
+    await ctx.report_progress(1, 2, "Fetching colormap...")
+    result = await client.explain_layer_colormap(layer_id)
+    await ctx.report_progress(2, 2, "Complete")
+    return result
 
 
 @mcp.tool(description="Fetch the pre-rendered legend image for a GIBS layer.")
 async def get_legend(
     layer_id: str,
+    ctx: Context,
     orientation: str = "horizontal",
 ):
     """Fetch the legend graphic for a layer.
@@ -494,14 +560,114 @@ async def get_legend(
         orientation: "horizontal" or "vertical".
     """
     client = _get_client()
+    await ctx.report_progress(1, 2, "Fetching legend...")
     legend_bytes = await client.get_legend_image(layer_id, orientation)
 
     if legend_bytes is None:
         return f"No legend available for '{layer_id}'."
 
+    await ctx.report_progress(2, 2, "Complete")
     return [f"Legend for {layer_id}", Image(data=legend_bytes, format="png")]
 
 
+@mcp.tool(
+    description="Get the exact data value at a geographic coordinate. "
+    "Reverse-maps the pixel color through the layer's colormap to recover "
+    "the numeric value with units. Only works for layers with colormaps."
+)
+async def query_point(
+    layer_id: str,
+    date: str,
+    lat: float,
+    lon: float,
+    ctx: Context,
+) -> str:
+    """Query the data value at a specific point.
+
+    Args:
+        layer_id: GIBS layer identifier (must have a colormap).
+        date: Date in YYYY-MM-DD format.
+        lat: Latitude in decimal degrees.
+        lon: Longitude in decimal degrees.
+    """
+    client = _get_client()
+
+    layer = client.get_layer(layer_id)
+    if layer is None:
+        return f"Layer '{layer_id}' not found."
+    if not layer.has_colormap:
+        return f"Layer '{layer_id}' has no colormap — cannot reverse-map pixel values."
+
+    await ctx.report_progress(1, 3, "Fetching pixel...")
+    result = await client.query_point(layer_id, date, lat, lon)
+    await ctx.report_progress(3, 3, "Complete")
+
+    return json.dumps(result, indent=2)
+
+
+@mcp.tool(
+    description="Fetch imagery for a layer across multiple dates. "
+    "Returns a sequence of images at regular intervals between start and end dates. "
+    "Useful for observing temporal progression (up to 12 frames)."
+)
+async def get_time_series(
+    layer_id: str,
+    start_date: str,
+    end_date: str,
+    ctx: Context,
+    bbox: list[float] | None = None,
+    place: str | None = None,
+    steps: int = 6,
+    width: int = 512,
+    height: int = 512,
+    projection: str = "auto",
+):
+    """Fetch imagery for multiple dates in a range.
+
+    Args:
+        layer_id: GIBS layer identifier.
+        start_date: First date (YYYY-MM-DD).
+        end_date: Last date (YYYY-MM-DD).
+        bbox: Bounding box as [west, south, east, north].
+        place: Place name to geocode.
+        steps: Number of evenly-spaced frames (2-12).
+        width: Image width in pixels.
+        height: Image height in pixels.
+        projection: EPSG code or "auto" to detect from bbox.
+    """
+    client = _get_client()
+
+    layer = client.get_layer(layer_id)
+    if layer is None:
+        return f"Layer '{layer_id}' not found."
+
+    steps = max(2, min(steps, 12))
+    dates = _generate_dates(start_date, end_date, steps)
+
+    await ctx.report_progress(1, steps + 2, "Resolving location...")
+    try:
+        resolved_bbox = await _resolve_bbox(client, bbox, place)
+    except Exception as exc:
+        return str(exc)
+
+    mime_type, short_fmt = _pick_format(layer)
+    epsg, proj_bbox = _pick_projection(resolved_bbox, projection)
+
+    results = []
+    for i, date in enumerate(dates):
+        await ctx.report_progress(
+            i + 2, steps + 2, f"Fetching {date} ({i + 1}/{len(dates)})..."
+        )
+        image_bytes = await client.get_wms_image(
+            layer_id, date, proj_bbox, width, height, mime_type, epsg
+        )
+        results.append(f"{layer.title} — {date}")
+        results.append(Image(data=image_bytes, format=short_fmt))
+
+    await ctx.report_progress(steps + 2, steps + 2, "Complete")
+    return results
+
+
 # ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
 #  TOOLS — Utility
 # ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
@@ -511,18 +677,20 @@ async def get_legend(
     description="Geocode a place name to geographic coordinates and bounding box. "
     "Uses OpenStreetMap Nominatim."
 )
-async def resolve_place(place: str) -> str:
+async def resolve_place(place: str, ctx: Context) -> str:
     """Resolve a place name to coordinates and bounding box.
 
     Args:
         place: Place name (e.g. "Tokyo", "California", "Great Barrier Reef").
     """
     client = _get_client()
+    await ctx.report_progress(1, 2, "Geocoding...")
     result = await client.resolve_place(place)
 
     if result is None:
         return f"Could not geocode '{place}'. Try a more specific name."
 
+    await ctx.report_progress(2, 2, "Complete")
     return json.dumps(
         {
             "display_name": result.display_name,
@@ -944,9 +1112,11 @@ def polar_watch(
     if pole.lower().startswith("ant"):
         region = "Antarctic"
         bbox = "[-180, -90, 180, -60]"
+        projection = "3031"
     else:
         region = "Arctic"
         bbox = "[-180, 60, 180, 90]"
+        projection = "3413"
 
     lines = [
         f"Monitor {region} ice and snow conditions on {date}.",
@@ -961,7 +1131,9 @@ def polar_watch(
         f"data exists for {date} on the selected layers.",
         "",
         f"3. **Polar imagery**: Fetch imagery using the {region} "
-        f"bounding box {bbox}. Get true color first, then ice/snow "
+        f"bounding box {bbox} with projection=\"{projection}\". "
+        f"This polar stereographic projection shows the region properly "
+        f"centered on the pole. Get true color first, then ice/snow "
         f"concentration layers.",
         "",
         "4. **Color interpretation**: Use explain_layer_colormap on each "
@@ -991,6 +1163,98 @@ def polar_watch(
     return "\n".join(lines)
 
 
+@mcp.prompt
+def quantitative_snapshot(
+    layer_id: str,
+    locations: str,
+    date: str,
+) -> str:
+    """Query exact data values at specific geographic coordinates.
+
+    Args:
+        layer_id: GIBS layer identifier with a colormap
+            (e.g. "AIRS_L3_Surface_Air_Temperature_Daily_Day").
+        locations: Comma-separated place names (e.g. "Tokyo, Sydney, Sahara Desert").
+        date: Date to query (YYYY-MM-DD).
+    """
+    location_list = [loc.strip() for loc in locations.split(",")]
+
+    lines = [
+        f"Query the exact {layer_id} values at these locations on {date}: "
+        f"{', '.join(location_list)}.",
+        "",
+        "Follow this workflow:",
+        "",
+        f"1. **Understand the scale**: Use explain_layer_colormap on "
+        f'"{layer_id}" to learn the value range, units, and what '
+        f"different colors represent. Fetch the legend with get_legend.",
+        "",
+        "2. **Resolve coordinates**: Use resolve_place for each location "
+        "to get precise lat/lon coordinates.",
+        "",
+        "3. **Query each point**: Use query_point for each location "
+        "with the resolved coordinates. This returns the exact data "
+        "value by reverse-mapping the pixel color through the colormap.",
+        "",
+        "4. **Compare and contextualize**: Present the results in a table "
+        "with location, value (in native and converted units if applicable), "
+        "and a qualitative description. Explain why the values differ — "
+        "latitude, altitude, season, land use, proximity to water, etc.",
+        "",
+        "5. **Verify with imagery**: Optionally fetch get_imagery for one "
+        "location to visually confirm the data pattern matches the "
+        "point query results.",
+    ]
+    return "\n".join(lines)
+
+
+@mcp.prompt
+def seasonal_timelapse(
+    layer_id: str,
+    location: str,
+    start_date: str,
+    end_date: str,
+) -> str:
+    """Watch a location change across multiple dates.
+
+    Args:
+        layer_id: GIBS layer identifier (e.g. "MODIS_Terra_NDVI_8Day").
+        location: Place to observe.
+        start_date: Beginning of the period (YYYY-MM-DD).
+        end_date: End of the period (YYYY-MM-DD).
+    """
+    lines = [
+        f"Show how {layer_id} changes at {location} "
+        f"from {start_date} to {end_date}.",
+        "",
+        "Follow this workflow:",
+        "",
+        f"1. **Layer info**: Use get_layer_info on \"{layer_id}\" to "
+        f"understand what it measures, its temporal period, and resolution.",
+        "",
+        "2. **Check dates**: Use check_layer_dates to verify data "
+        f"covers {start_date} through {end_date}.",
+        "",
+        "3. **Color interpretation**: Use explain_layer_colormap to "
+        "understand the value encoding, so you can describe changes "
+        "in scientific terms rather than just color names.",
+        "",
+        f"4. **Time series**: Use get_time_series to fetch imagery "
+        f'for "{location}" from {start_date} to {end_date}. '
+        f"This returns a sequence of images at regular intervals.",
+        "",
+        "5. **Narrate the progression**: For each frame, describe "
+        "what's visible and how it differs from the previous frame. "
+        "Reference colormap values to quantify the change. "
+        "Identify the overall trend and any notable transitions.",
+        "",
+        "6. **Bookend comparison**: Use compare_dates with "
+        f"{start_date} and {end_date} for a direct side-by-side "
+        "showing the total change across the period.",
+    ]
+    return "\n".join(lines)
+
+
 # ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
 #  Entry point
 # ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━

tests/test_geo.py

@@ -3,7 +3,7 @@
 import httpx
 import respx
 
-from mcgibs.geo import bbox_from_point, expand_bbox, geocode
+from mcgibs.geo import bbox_from_point, bbox_to_polar, expand_bbox, geocode
 from mcgibs.models import BBox
 
 # ---------------------------------------------------------------------------
@@ -149,3 +149,45 @@ def test_bbox_from_point_clamping():
     assert bbox.south == 87.5
     assert bbox.west == -2.0
     assert bbox.east == 2.0
+
+
+# ---------------------------------------------------------------------------
+# bbox_to_polar() tests
+# ---------------------------------------------------------------------------
+
+
+def test_polar_full_arctic_extent():
+    """Full Arctic bbox should return the GIBS documented extent."""
+    bbox = BBox(west=-180, south=60, east=180, north=90)
+    xmin, ymin, xmax, ymax = bbox_to_polar(bbox, "3413")
+
+    assert xmin == -4194304.0
+    assert ymin == -4194304.0
+    assert xmax == 4194304.0
+    assert ymax == 4194304.0
+
+
+def test_polar_full_antarctic_extent():
+    """Full Antarctic bbox should return the GIBS documented extent."""
+    bbox = BBox(west=-180, south=-90, east=180, north=-60)
+    xmin, _ymin, xmax, _ymax = bbox_to_polar(bbox, "3031")
+
+    assert xmin == -4194304.0
+    assert xmax == 4194304.0
+
+
+def test_polar_subregion_is_bounded():
+    """A small Arctic subregion should produce coordinates smaller than full extent."""
+    bbox = BBox(west=-60, south=70, east=-20, north=80)
+    xmin, ymin, xmax, ymax = bbox_to_polar(bbox, "3413")
+
+    assert -4194304 < xmin < xmax < 4194304
+    assert -4194304 < ymin < ymax < 4194304
+
+
+def test_polar_passthrough_for_4326():
+    """Non-polar EPSG should return geographic coords unchanged."""
+    bbox = BBox(west=-120, south=30, east=-110, north=40)
+    result = bbox_to_polar(bbox, "4326")
+
+    assert result == (-120, 30, -110, 40)

tests/test_tools.py

@@ -166,6 +166,8 @@ async def test_list_tools(capabilities_xml):
                 "get_legend",
                 "resolve_place",
                 "build_tile_url",
+                "query_point",
+                "get_time_series",
             }
 
             for name in expected:
@@ -303,6 +305,8 @@ async def test_list_prompts(capabilities_xml):
                 "layer_deep_dive",
                 "multi_layer_story",
                 "polar_watch",
+                "quantitative_snapshot",
+                "seasonal_timelapse",
             }
 
             for name in expected:
@@ -358,3 +362,191 @@ async def test_colormap_resource(capabilities_xml, colormap_xml):
     finally:
         await server_module._client.close()
         server_module._client = None
+
+
+# ---------------------------------------------------------------------------
+# query_point tests
+# ---------------------------------------------------------------------------
+
+
+def _make_fake_png(rgb: tuple[int, int, int], size: int = 3) -> bytes:
+    """Create a tiny PNG with uniform color for mocking point queries."""
+    buf = BytesIO()
+    PILImage.new("RGBA", (size, size), (*rgb, 255)).save(buf, format="PNG")
+    return buf.getvalue()
+
+
+@respx.mock
+async def test_query_point(capabilities_xml, colormap_xml):
+    """query_point reverse-maps pixel RGB through colormap to data value."""
+    respx.get(url__regex=r".*WMTSCapabilities\.xml").mock(
+        return_value=httpx.Response(200, text=capabilities_xml)
+    )
+    # Return a 3x3 PNG with rgb (255,100,50) → maps to [290,300) K
+    respx.get(url__regex=r".*wms\.cgi.*").mock(
+        return_value=httpx.Response(
+            200,
+            content=_make_fake_png((255, 100, 50)),
+            headers={"content-type": "image/png"},
+        )
+    )
+    respx.get(url__regex=r".*colormaps/v1\.3/.*\.xml").mock(
+        return_value=httpx.Response(200, text=colormap_xml)
+    )
+
+    server_module._client = await _init_mock_client(capabilities_xml)
+    try:
+        async with Client(mcp) as client:
+            result = await client.call_tool(
+                "query_point",
+                {
+                    "layer_id": "AIRS_L3_Surface_Air_Temperature_Daily_Day",
+                    "date": "2025-06-01",
+                    "lat": 35.67,
+                    "lon": 139.65,
+                },
+            )
+            data = json.loads(result.content[0].text)
+
+            assert data["lat"] == 35.67
+            assert data["lon"] == 139.65
+            assert data["quality"] == "exact"  # exact RGB match
+            assert "value" in data
+    finally:
+        await server_module._client.close()
+        server_module._client = None
+
+
+@respx.mock
+async def test_query_point_no_colormap(capabilities_xml):
+    """query_point returns error for layers without colormaps."""
+    respx.get(url__regex=r".*WMTSCapabilities\.xml").mock(
+        return_value=httpx.Response(200, text=capabilities_xml)
+    )
+
+    server_module._client = await _init_mock_client(capabilities_xml)
+    try:
+        async with Client(mcp) as client:
+            result = await client.call_tool(
+                "query_point",
+                {
+                    "layer_id": "MODIS_Terra_CorrectedReflectance_TrueColor",
+                    "date": "2025-06-01",
+                    "lat": 35.67,
+                    "lon": 139.65,
+                },
+            )
+            text = result.content[0].text
+            assert "no colormap" in text
+    finally:
+        await server_module._client.close()
+        server_module._client = None
+
+
+# ---------------------------------------------------------------------------
+# get_time_series tests
+# ---------------------------------------------------------------------------
+
+
+@respx.mock
+async def test_get_time_series(capabilities_xml):
+    """get_time_series returns multiple dated images."""
+    respx.get(url__regex=r".*WMTSCapabilities\.xml").mock(
+        return_value=httpx.Response(200, text=capabilities_xml)
+    )
+    respx.get(url__regex=r".*wms\.cgi.*").mock(
+        return_value=httpx.Response(
+            200,
+            content=_make_fake_jpeg(),
+            headers={"content-type": "image/jpeg"},
+        )
+    )
+
+    server_module._client = await _init_mock_client(capabilities_xml)
+    try:
+        async with Client(mcp) as client:
+            result = await client.call_tool(
+                "get_time_series",
+                {
+                    "layer_id": "MODIS_Terra_CorrectedReflectance_TrueColor",
+                    "start_date": "2025-01-01",
+                    "end_date": "2025-06-01",
+                    "bbox": [-120.0, 30.0, -110.0, 40.0],
+                    "steps": 3,
+                },
+            )
+            # 3 steps → 3 text labels + 3 images = 6 content items
+            texts = [c for c in result.content if c.type == "text"]
+            images = [c for c in result.content if c.type == "image"]
+
+            assert len(texts) == 3
+            assert len(images) == 3
+            assert "2025-01-01" in texts[0].text
+            assert "2025-06-01" in texts[-1].text
+    finally:
+        await server_module._client.close()
+        server_module._client = None
+
+
+# ---------------------------------------------------------------------------
+# New prompt tests
+# ---------------------------------------------------------------------------
+
+
+@respx.mock
+async def test_quantitative_snapshot_prompt(capabilities_xml):
+    """quantitative_snapshot prompt includes locations and tool references."""
+    respx.get(url__regex=r".*WMTSCapabilities\.xml").mock(
+        return_value=httpx.Response(200, text=capabilities_xml)
+    )
+
+    server_module._client = await _init_mock_client(capabilities_xml)
+    try:
+        async with Client(mcp) as client:
+            result = await client.get_prompt(
+                "quantitative_snapshot",
+                {
+                    "layer_id": "AIRS_L3_Surface_Air_Temperature_Daily_Day",
+                    "locations": "Tokyo, Sydney",
+                    "date": "2025-06-01",
+                },
+            )
+            text = result.messages[0].content.text
+
+            assert "Tokyo" in text
+            assert "Sydney" in text
+            assert "query_point" in text
+            assert "explain_layer_colormap" in text
+    finally:
+        await server_module._client.close()
+        server_module._client = None
+
+
+@respx.mock
+async def test_seasonal_timelapse_prompt(capabilities_xml):
+    """seasonal_timelapse prompt includes dates and get_time_series reference."""
+    respx.get(url__regex=r".*WMTSCapabilities\.xml").mock(
+        return_value=httpx.Response(200, text=capabilities_xml)
+    )
+
+    server_module._client = await _init_mock_client(capabilities_xml)
+    try:
+        async with Client(mcp) as client:
+            result = await client.get_prompt(
+                "seasonal_timelapse",
+                {
+                    "layer_id": "MODIS_Terra_CorrectedReflectance_TrueColor",
+                    "location": "Amazon Rainforest",
+                    "start_date": "2025-01-01",
+                    "end_date": "2025-12-01",
+                },
+            )
+            text = result.messages[0].content.text
+
+            assert "Amazon Rainforest" in text
+            assert "get_time_series" in text
+            assert "2025-01-01" in text
+            assert "2025-12-01" in text
+    finally:
+        await server_module._client.close()
+        server_module._client = None

uv.lock

@@ -594,7 +594,7 @@ wheels = [
 
 [[package]]
 name = "mcgibs"
-version = "2026.2.18"
+version = "2026.2.19"
 source = { editable = "." }
 dependencies = [
     { name = "defusedxml" },