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],
|
|||
|
|
}
|