99 lines
3.5 KiB
Python
99 lines
3.5 KiB
Python
"""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],
|
||
}
|