Compare commits

..

No commits in common. "main" and "v2026.02.19" have entirely different histories.

12 changed files with 46 additions and 247 deletions

View File

@ -311,13 +311,6 @@ The assistant fetches a sequence of images at regular intervals:
[6 images: Mar 1, Apr 13, May 25, Jul 6, Aug 17, Sep 1] [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: > **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. > 1. **March** — The Plains are mostly tan/brown (NDVI ~0.1), winter dormancy. Only southern Texas shows hints of green.
@ -338,7 +331,7 @@ The assistant fetches a sequence of images at regular intervals:
The `polar_watch` prompt guides the assistant to use polar stereographic projection: The `polar_watch` prompt guides the assistant to use polar stereographic projection:
``` ```
→ get_imagery("AMSRU2_Sea_Ice_Concentration_12km", "2025-03-01", → get_imagery("AMSRU2_Sea_Ice_Concentration_12km", "2025-09-01",
bbox=[-180, 60, 180, 90], bbox=[-180, 60, 180, 90],
projection="3413") projection="3413")
[image: Arctic sea ice in polar stereographic — circular view centered on the pole] [image: Arctic sea ice in polar stereographic — circular view centered on the pole]
@ -350,17 +343,12 @@ The `polar_watch` prompt guides the assistant to use polar stereographic project
- light grey = ~70% - light grey = ~70%
- white = ~100% (solid pack ice) - white = ~100% (solid pack ice)
→ build_tile_url("AMSRU2_Sea_Ice_Concentration_12km", "2025-03-01", → build_tile_url("AMSRU2_Sea_Ice_Concentration_12km", "2025-09-01",
zoom=2, row=0, col=0, projection="3413") zoom=2, row=0, col=0, projection="3413")
https://gibs.earthdata.nasa.gov/wmts/epsg3413/best/AMSRU2_Sea_Ice_Concentration_12km/... https://gibs.earthdata.nasa.gov/wmts/epsg3413/best/AMSRU2_Sea_Ice_Concentration_12km/...
``` ```
<p align="center"> > **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.
<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 ## Projections

Binary file not shown.

Before

Width:  |  Height:  |  Size: 118 KiB

After

Width:  |  Height:  |  Size: 190 KiB

Binary file not shown.

Before

Width:  |  Height:  |  Size: 80 KiB

Binary file not shown.

Before

Width:  |  Height:  |  Size: 663 KiB

Binary file not shown.

Before

Width:  |  Height:  |  Size: 607 KiB

Binary file not shown.

Before

Width:  |  Height:  |  Size: 686 KiB

View File

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

View File

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

View File

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

View File

@ -1,7 +1,6 @@
"""Async Nominatim geocoding with bbox utilities.""" """Async Nominatim geocoding with bbox utilities."""
import logging import logging
import math
import httpx import httpx
@ -130,99 +129,3 @@ def bbox_from_point(lat: float, lon: float, radius_deg: float = 0.5) -> BBox:
east=min(180.0, lon + radius_deg), east=min(180.0, lon + radius_deg),
north=min(90.0, lat + 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))

View File

@ -23,7 +23,7 @@ from fastmcp.utilities.types import Image
from mcgibs.capabilities import search_layers from mcgibs.capabilities import search_layers
from mcgibs.client import GIBSClient from mcgibs.client import GIBSClient
from mcgibs.constants import PROJECTIONS from mcgibs.constants import PROJECTIONS
from mcgibs.geo import bbox_to_polar, expand_bbox from mcgibs.geo import expand_bbox
from mcgibs.models import BBox from mcgibs.models import BBox
log = logging.getLogger(__name__) log = logging.getLogger(__name__)
@ -304,43 +304,15 @@ async def _resolve_bbox(
raise ValueError("Provide either 'bbox' [west, south, east, north] or 'place' name.") raise ValueError("Provide either 'bbox' [west, south, east, north] or 'place' name.")
def _pick_projection(bbox: BBox, projection: str) -> tuple[str, BBox]: def _pick_projection(bbox: BBox, projection: str) -> str:
"""Resolve projection and transform bbox if needed. """Resolve projection — auto-detect polar if set to 'auto'."""
if projection != "auto":
Returns (epsg, transformed_bbox). For polar projections the bbox is return projection
converted from geographic coordinates to projected meters.
"""
if projection == "auto":
if bbox.south >= 55: if bbox.south >= 55:
projection = "3413" return "3413"
elif bbox.north <= -55: if bbox.north <= -55:
projection = "3031" return "3031"
else: return "4326"
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]: def _generate_dates(start: str, end: str, steps: int) -> list[str]:
@ -367,7 +339,7 @@ async def get_imagery(
place: str | None = None, place: str | None = None,
width: int = 1024, width: int = 1024,
height: int = 1024, height: int = 1024,
format: str = "auto", format: str = "jpeg",
projection: str = "auto", projection: str = "auto",
): ):
"""Fetch GIBS imagery via WMS. """Fetch GIBS imagery via WMS.
@ -379,7 +351,7 @@ async def get_imagery(
place: Place name to geocode (e.g. "California", "Amazon Basin"). place: Place name to geocode (e.g. "California", "Amazon Basin").
width: Image width in pixels. width: Image width in pixels.
height: Image height in pixels. height: Image height in pixels.
format: Image format "jpeg", "png", or "auto" (picks PNG for colormapped layers). format: Image format "jpeg" or "png".
projection: EPSG code ("4326", "3857", "3413", "3031") or "auto" to detect from bbox. projection: EPSG code ("4326", "3857", "3413", "3031") or "auto" to detect from bbox.
""" """
client = _get_client() client = _get_client()
@ -388,22 +360,23 @@ async def get_imagery(
if layer is None: if layer is None:
return f"Layer '{layer_id}' not found." 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...") await ctx.report_progress(1, 4, "Resolving location...")
try: try:
resolved_bbox = await _resolve_bbox(client, bbox, place) resolved_bbox = await _resolve_bbox(client, bbox, place)
except Exception as exc: except Exception as exc:
return str(exc) return str(exc)
epsg, proj_bbox = _pick_projection(resolved_bbox, projection) epsg = _pick_projection(resolved_bbox, projection)
await ctx.report_progress(2, 4, "Fetching imagery from GIBS...") await ctx.report_progress(2, 4, "Fetching imagery from GIBS...")
image_bytes = await client.get_wms_image( image_bytes = await client.get_wms_image(
layer_id, date, proj_bbox, width, height, mime_type, epsg, layer_id,
date,
resolved_bbox,
width,
height,
f"image/{format}",
epsg,
) )
description = ( description = (
@ -415,7 +388,7 @@ async def get_imagery(
description += "\nTip: use explain_layer_colormap to understand what the colors represent." description += "\nTip: use explain_layer_colormap to understand what the colors represent."
await ctx.report_progress(4, 4, "Complete") await ctx.report_progress(4, 4, "Complete")
return [description, Image(data=image_bytes, format=short_fmt)] return [description, Image(data=image_bytes, format=format)]
@mcp.tool( @mcp.tool(
@ -453,14 +426,17 @@ async def compare_dates(
except Exception as exc: except Exception as exc:
return str(exc) return str(exc)
mime_type, short_fmt = _pick_format(layer) epsg = _pick_projection(resolved_bbox, projection)
epsg, proj_bbox = _pick_projection(resolved_bbox, projection)
await ctx.report_progress(2, 5, f"Fetching imagery for {date_before}...") await ctx.report_progress(2, 5, f"Fetching imagery for {date_before}...")
await ctx.report_progress(3, 5, f"Fetching imagery for {date_after}...") await ctx.report_progress(3, 5, f"Fetching imagery for {date_after}...")
composite_bytes = await client.compare_dates( composite_bytes = await client.compare_dates(
layer_id, date_before, date_after, proj_bbox, layer_id,
image_format=mime_type, epsg=epsg, date_before,
date_after,
resolved_bbox,
image_format="image/jpeg",
epsg=epsg,
) )
description = ( description = (
@ -470,7 +446,7 @@ async def compare_dates(
) )
await ctx.report_progress(5, 5, "Complete") await ctx.report_progress(5, 5, "Complete")
return [description, Image(data=composite_bytes, format=short_fmt)] return [description, Image(data=composite_bytes, format="jpeg")]
@mcp.tool( @mcp.tool(
@ -509,12 +485,17 @@ async def get_imagery_composite(
except Exception as exc: except Exception as exc:
return str(exc) return str(exc)
epsg, proj_bbox = _pick_projection(resolved_bbox, projection) epsg = _pick_projection(resolved_bbox, projection)
layer_names = ", ".join(layer_ids) layer_names = ", ".join(layer_ids)
await ctx.report_progress(2, 3, f"Fetching {len(layer_ids)}-layer composite...") await ctx.report_progress(2, 3, f"Fetching {len(layer_ids)}-layer composite...")
image_bytes = await client.get_wms_composite( image_bytes = await client.get_wms_composite(
layer_ids, date, proj_bbox, width, height, epsg=epsg, layer_ids,
date,
resolved_bbox,
width,
height,
epsg=epsg,
) )
description = ( description = (
@ -650,8 +631,7 @@ async def get_time_series(
except Exception as exc: except Exception as exc:
return str(exc) return str(exc)
mime_type, short_fmt = _pick_format(layer) epsg = _pick_projection(resolved_bbox, projection)
epsg, proj_bbox = _pick_projection(resolved_bbox, projection)
results = [] results = []
for i, date in enumerate(dates): for i, date in enumerate(dates):
@ -659,10 +639,10 @@ async def get_time_series(
i + 2, steps + 2, f"Fetching {date} ({i + 1}/{len(dates)})..." i + 2, steps + 2, f"Fetching {date} ({i + 1}/{len(dates)})..."
) )
image_bytes = await client.get_wms_image( image_bytes = await client.get_wms_image(
layer_id, date, proj_bbox, width, height, mime_type, epsg layer_id, date, resolved_bbox, width, height, "image/jpeg", epsg
) )
results.append(f"{layer.title}{date}") results.append(f"{layer.title}{date}")
results.append(Image(data=image_bytes, format=short_fmt)) results.append(Image(data=image_bytes, format="jpeg"))
await ctx.report_progress(steps + 2, steps + 2, "Complete") await ctx.report_progress(steps + 2, steps + 2, "Complete")
return results return results

View File

@ -3,7 +3,7 @@
import httpx import httpx
import respx import respx
from mcgibs.geo import bbox_from_point, bbox_to_polar, expand_bbox, geocode from mcgibs.geo import bbox_from_point, expand_bbox, geocode
from mcgibs.models import BBox from mcgibs.models import BBox
# --------------------------------------------------------------------------- # ---------------------------------------------------------------------------
@ -149,45 +149,3 @@ def test_bbox_from_point_clamping():
assert bbox.south == 87.5 assert bbox.south == 87.5
assert bbox.west == -2.0 assert bbox.west == -2.0
assert bbox.east == 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)