添加火星TMS服务

app_mars.py

@@ -0,0 +1,98 @@
+"""rio-tiler dynamic tile server for Mars MGS MOLA ClrShade global basemap"""
+
+import math
+from pathlib import Path
+
+from fastapi import FastAPI, HTTPException
+from starlette.requests import Request
+from starlette.responses import Response
+
+from pyproj import CRS
+from morecantile import TileMatrixSet
+
+from rio_tiler.errors import TileOutsideBounds
+from rio_tiler.io import Reader
+from rio_tiler.profiles import img_profiles
+
+DATA_PATH = str(
+    Path(__file__).parent / "data" / "Mars_MGS_MOLA_ClrShade_merge_global_463m.tif"
+)
+
+# TIF 的原生 CRS 是 SimpleCylindrical Mars（等距柱状投影，单位：米）
+# 用 proj4 字符串定义，直接匹配 TIF —— 参考 lunaserv-python 的 proj4 字符串方式
+#
+# 关键设计决策：TMS 与 TIF 使用同一个 CRS（eqc 米制），rio-tiler 读取时无需任何跨 CRS 转换。
+# 若改用 longlat（度），PROJ 会尝试通过 WGS84 做 datum shift：火星球体(R=3396km) vs
+# 地球椭球(R=6378km) 差异极大，导致中间计算出现 "PROJ: eqc: Invalid latitude" 错误。
+MARS_CRS = CRS.from_proj4(
+    "+proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +a=3396190 +b=3396190 +units=m +no_defs"
+)
+
+# 火星球体在等距柱状投影下的全球范围（米）
+# 等距柱状: x = R × lon_rad, y = R × lat_rad
+# → max_x = R × π（对应 180°），max_y = R × π/2（对应 90°）
+_R = 3396190
+_MAX_X = _R * math.pi        # ≈ 10,669,442 m
+_MAX_Y = _R * math.pi / 2    # ≈  5,334,721 m
+
+# WorldCRS84Quad 等效 TMS（米制版本）：
+# - 范围覆盖全球等距柱状投影
+# - matrix_scale=[2, 1] → zoom 0 有 2×1 个瓦片，与 Cesium GeographicTilingScheme(2×1) 对齐
+# - 瓦片编号 (z, x, y) 与地理位置的对应关系和经纬度版完全一致，Cesium 无需任何修改
+MARS_TMS = TileMatrixSet.custom(
+    crs=MARS_CRS,
+    extent=(-_MAX_X, -_MAX_Y, _MAX_X, _MAX_Y),
+    matrix_scale=[2, 1],
+)
+
+app = FastAPI(
+    title="rio-tiler Mars",
+    description="Mars MGS MOLA Color-Shaded Relief global tile server",
+)
+
+
+@app.get(
+    "/tiles/{z}/{x}/{y}.png",
+    responses={
+        200: {"content": {"image/png": {}}, "description": "Return a PNG tile."},
+        404: {"description": "Tile outside bounds."},
+    },
+    response_class=Response,
+    description="Read Mars tile and return PNG (uint8 RGB, no rescaling needed)",
+)
+def tile(z: int, x: int, y: int):
+    """Return a map tile for the given z/x/y coordinates."""
+    try:
+        with Reader(DATA_PATH, tms=MARS_TMS) as cog:
+            img = cog.tile(x, y, z)
+    except TileOutsideBounds:
+        raise HTTPException(status_code=404, detail="Tile outside bounds")
+
+    content = img.render(
+        img_format="PNG",
+        **img_profiles.get("png", {}),
+    )
+    return Response(content, media_type="image/png")
+
+
+@app.get("/tilejson.json", responses={200: {"description": "Return a TileJSON document"}})
+def tilejson(request: Request):
+    """Return a TileJSON 2.2.0 document."""
+    with Reader(DATA_PATH, tms=MARS_TMS) as cog:
+        bounds = cog.get_geographic_bounds(cog.tms.rasterio_geographic_crs)
+        minzoom = cog.minzoom
+        maxzoom = cog.maxzoom
+
+    base_url = str(request.base_url).rstrip("/")
+    tile_url = f"{base_url}/tiles/{{z}}/{{x}}/{{y}}.png"
+
+    return {
+        "tilejson": "2.2.0",
+        "name": "Mars MGS MOLA ClrShade",
+        "description": "Mars MGS MOLA Color-Shaded Relief Global Mosaic 463m/px",
+        "tiles": [tile_url],
+        "bounds": list(bounds),
+        "minzoom": minzoom,
+        "maxzoom": maxzoom,
+        "center": [0, 0, minzoom],
+    }