stac_mosaic_api.py

@@ -67,12 +67,10 @@ def bbox_mosaic(query: BBoxQuery):
         if not image_paths:
             raise HTTPException(status_code=404, detail="未查询到图像")
 
-        # TODO 这个函数的输入可以用*args, **kwargs，指定了max_size，后面记得改
         def part_reader(src_path, part):
             with Reader(src_path) as cog:
                 return cog.part(part, max_size=1024)
 
-        # TODO 目前部分读取和镶嵌一起做做的，以后做实验需要解耦
         img_buf = render_mosaic(image_paths, part_reader, query.bbox)
         return StreamingResponse(img_buf, media_type="image/png")
     except Exception as e:
@@ -90,8 +88,7 @@ def polygon_mosaic(polygon: GeoJSONPolygon = Body(...)):
             with Reader(src_path) as cog:
                 return cog.feature(feat, max_size=1024)
 
-        img_buf = render_mosaic(
+        img_buf = render_mosaic(image_paths, feature_reader, polygon.model_dump())
-            image_paths, feature_reader, polygon.model_dump())
         return StreamingResponse(img_buf, media_type="image/png")
     except Exception as e:
         raise HTTPException(status_code=500, detail=str(e))

stac_mosaic_stream_api.py

@@ -1,114 +0,0 @@
-from fastapi import FastAPI, Query, Body, HTTPException
-from typing import List
-from pydantic import BaseModel
-from shapely.geometry import Polygon, mapping, shape
-from rio_tiler.io import Reader
-from rio_tiler.mosaic import mosaic_reader
-import requests
-from fastapi.responses import StreamingResponse
-from io import BytesIO
-import morecantile
-from starlette.responses import StreamingResponse
-
-# FastAPI 实例
-app = FastAPI(title="STAC Mosaic API")
-
-# STAC 配置
-STAC_API_URL = "http://localhost:8082"
-COLLECTION_ID = "geosat1"
-
-
-# ---------- 数据模型 ----------
-class BBoxQuery(BaseModel):
-    bbox: List[float]  # [minx, miny, maxx, maxy]
-
-
-class GeoJSONPolygon(BaseModel):
-    type: str
-    coordinates: List
-
-
-# ---------- 实用函数 ----------
-def fetch_tiles_by_polygon(polygon: GeoJSONPolygon) -> List[tuple]:
-    """
-    根据多边形获取相交的瓦片列表
-    """
-    zoom = 11
-    tms = morecantile.tms.get("WebMercatorQuad")
-    # 将 GeoJSONPolygon 转为 shapely Polygon
-    shapely_poly = shape(polygon.model_dump())
-    tiles = list(tms.tiles(*shapely_poly.bounds, zoom))
-
-    def tile_polygon(x, y, z):
-        bounds = tms.bounds(x, y, z)
-        return Polygon([
-            (bounds.left, bounds.bottom),
-            (bounds.left, bounds.top),
-            (bounds.right, bounds.top),
-            (bounds.right, bounds.bottom),
-            (bounds.left, bounds.bottom)
-        ])
-
-    intersecting_tiles = []
-    for x, y, z in tiles:
-        if shapely_poly.intersects(tile_polygon(x, y, z)):
-            intersecting_tiles.append((x, y, z))
-
-    return intersecting_tiles
-
-
-def fetch_items_by_tile(x, y, z) -> List[str]:
-    """
-    根据瓦片坐标 (x, y, z) 查询与该瓦片 bbox 相交的影像
-    """
-    tms = morecantile.tms.get("WebMercatorQuad")
-    bounds = tms.bounds(x, y, z)
-    bbox = [bounds.left, bounds.bottom, bounds.right, bounds.top]
-    url = f"{STAC_API_URL}/collections/{COLLECTION_ID}/items"
-    params = {"bbox": ",".join(map(str, bbox))}
-    r = requests.get(url, params=params)
-    r.raise_for_status()
-    data = r.json()
-    return [feat["assets"]["image"]["href"] for feat in data.get("features", [])]
-
-
-def tile_mosaic_buf(x, y, z):
-    image_paths = fetch_items_by_tile(x, y, z)
-    if not image_paths:
-        raise HTTPException(status_code=404, detail="未查询到图像")
-
-    def tile_reader(src_path, x, y, z):
-        with Reader(src_path) as cog:
-            return cog.tile(x, y, z, tilesize=256)
-    img, _ = mosaic_reader(image_paths, tile_reader, x, y, z)
-    buf = BytesIO(img.render(img_format="PNG"))
-    buf.seek(0)
-    return buf
-
-
-# ---------- API 路由 ----------
-@app.post("/polygon_tiles_mosaic_stream", summary="逐个返回瓦片镶嵌结果")
-def polygon_tiles_mosaic_stream(polygon: GeoJSONPolygon = Body(...)):
-    def tile_generator():
-        intersecting_tiles = fetch_tiles_by_polygon(polygon)
-        for x, y, z in intersecting_tiles:
-            buf = tile_mosaic_buf(x, y, z)
-            if buf:
-                # 每个 tile 之间用自定义分隔符分割，前端需配合解析
-                yield b"--tile-boundary--\n"
-                yield buf.read()
-
-            # TODO 下面的代码是可选的，如果需要返回瓦片的元数据，可以取消注释
-            # if buf:
-            #     img_bytes = buf.getvalue()
-            #     meta = {
-            #         "x": x,
-            #         "y": y,
-            #         "z": z,
-            #         "content_type": "image/png",
-            #         "length": len(img_bytes)
-            #     }
-            #     yield b"--tile-boundary--\n"
-            #     yield (json.dumps(meta) + "\n").encode("utf-8")
-            #     yield img_bytes
-    return StreamingResponse(tile_generator(), media_type="application/octet-stream")

tools/viusal_tifs_gps.py

@@ -1,39 +0,0 @@
-import rasterio
-from rasterio.warp import transform
-import matplotlib.pyplot as plt
-import glob
-import os
-
-def get_bounds_lonlat(dataset):
-    # 直接用元数据获取边界
-    left, bottom, right, top = dataset.bounds
-    # 四角顺序：左下、左上、右上、右下
-    xs = [left, left, right, right]
-    ys = [bottom, top, top, bottom]
-    # 转为经纬度
-    lons, lats = transform(dataset.crs, 'EPSG:4326', xs, ys)
-    return list(zip(lons, lats))
-
-if __name__ == "__main__":
-    tif_dir = "test_geosat1_imgs"
-    pattern = os.path.join(tif_dir, "Jilin*.tif")
-    all_bounds = []
-    for tif in glob.glob(pattern):
-        with rasterio.open(tif) as ds:
-            corners_lonlat = get_bounds_lonlat(ds)
-            all_bounds.append((os.path.basename(tif), corners_lonlat))
-
-    plt.figure(figsize=(10, 10))
-    for name, corners in all_bounds:
-        # 闭合多边形
-        corners_closed = corners + [corners[0]]
-        lons, lats = zip(*corners_closed)
-        plt.plot(lons, lats, marker='o', label=name)
-        for lon, lat in corners:
-            plt.text(lon, lat, f'({lon:.4f},{lat:.4f})', color='red', fontsize=8)
-    plt.xlabel("Longitude")
-    plt.ylabel("Latitude")
-    plt.title("Jilin* 影像边界经纬度分布")
-    plt.legend()
-    plt.grid(True)
-    plt.show()