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

@@ -311,6 +311,13 @@ The assistant fetches a sequence of images at regular intervals:
   [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.
@@ -331,7 +338,7 @@ The assistant fetches a sequence of images at regular intervals:
 The `polar_watch` prompt guides the assistant to use polar stereographic projection:
 
 ```
-→ get_imagery("AMSRU2_Sea_Ice_Concentration_12km", "2025-09-01",
+→ 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]
@@ -343,12 +350,17 @@ The `polar_watch` prompt guides the assistant to use polar stereographic project
     - light grey = ~70%
     - white = ~100% (solid pack ice)
 
-→ build_tile_url("AMSRU2_Sea_Ice_Concentration_12km", "2025-09-01",
+→ 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/...
 ```
 
-> **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 September minimum is clearly visible: solid pack ice (white, 100%) is confined to a shrinking core, surrounded by a wide margin of partial ice (blues and cyans, 15-40%). The embeddable WMTS tile URL is ready for dashboards or reports.
+<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
 

pyproject.toml

@@ -1,6 +1,6 @@
 [project]
 name = "mcgibs"
-version = "2026.02.19"
+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

@@ -98,7 +98,37 @@ async def main():
             )
             explanation = result.content[0].text
             (OUT / "temperature-colormap.txt").write_text(explanation)
-            print(f"  Saved colormap explanation", file=sys.stderr)
+            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()

src/mcgibs/__init__.py

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

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

@@ -23,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__)
@@ -304,15 +304,43 @@ async def _resolve_bbox(
     raise ValueError("Provide either 'bbox' [west, south, east, north] or 'place' name.")
 
 
-def _pick_projection(bbox: BBox, projection: str) -> str:
+def _pick_projection(bbox: BBox, projection: str) -> tuple[str, BBox]:
-    """Resolve projection — auto-detect polar if set to 'auto'."""
+    """Resolve projection and transform bbox if needed.
-    if projection != "auto":
+
-        return projection
+    Returns (epsg, transformed_bbox). For polar projections the bbox is
+    converted from geographic coordinates to projected meters.
+    """
+    if projection == "auto":
         if bbox.south >= 55:
-        return "3413"
+            projection = "3413"
-    if bbox.north <= -55:
+        elif bbox.north <= -55:
-        return "3031"
+            projection = "3031"
-    return "4326"
+        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]:
@@ -339,7 +367,7 @@ 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.
@@ -351,7 +379,7 @@ 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()
@@ -360,23 +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 = _pick_projection(resolved_bbox, projection)
+    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}",
-        epsg,
     )
 
     description = (
@@ -388,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(
@@ -426,17 +453,14 @@ async def compare_dates(
     except Exception as exc:
         return str(exc)
 
-    epsg = _pick_projection(resolved_bbox, projection)
+    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,
-        image_format="image/jpeg",
-        epsg=epsg,
     )
 
     description = (
@@ -446,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(
@@ -485,17 +509,12 @@ async def get_imagery_composite(
     except Exception as exc:
         return str(exc)
 
-    epsg = _pick_projection(resolved_bbox, projection)
+    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,
-        epsg=epsg,
     )
 
     description = (
@@ -631,7 +650,8 @@ async def get_time_series(
     except Exception as exc:
         return str(exc)
 
-    epsg = _pick_projection(resolved_bbox, projection)
+    mime_type, short_fmt = _pick_format(layer)
+    epsg, proj_bbox = _pick_projection(resolved_bbox, projection)
 
     results = []
     for i, date in enumerate(dates):
@@ -639,10 +659,10 @@ async def get_time_series(
             i + 2, steps + 2, f"Fetching {date} ({i + 1}/{len(dates)})..."
         )
         image_bytes = await client.get_wms_image(
-            layer_id, date, resolved_bbox, width, height, "image/jpeg", epsg
+            layer_id, date, proj_bbox, width, height, mime_type, epsg
         )
         results.append(f"{layer.title} — {date}")
-        results.append(Image(data=image_bytes, format="jpeg"))
+        results.append(Image(data=image_bytes, format=short_fmt))
 
     await ctx.report_progress(steps + 2, steps + 2, "Complete")
     return results

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)