python-fastapi-mm-zhcs-yj-api/routers/api/spatialAnalysis/point.py

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

from fastapi import APIRouter, Request, Depends, Query, HTTPException, status,Path
from common.security import valid_access_token
from fastapi.responses import JSONResponse,Response
from fastapi.responses import StreamingResponse
from sqlalchemy.orm import Session
from sqlalchemy import and_, or_,text,literal
from sqlalchemy.sql import func
from sqlalchemy.future import select
from common.auth_user import *
from pydantic import BaseModel
from typing import Any, Dict
# import contextily as ctx
# import geopandas as gpd
# from matplotlib import pyplot as plt
import io
from database import get_db
from typing import List
from models import *
from utils import *
from utils.ry_system_util import *
from utils.video_util import *
from collections import defaultdict
import traceback
from concurrent.futures import ThreadPoolExecutor, as_completed
from multiprocessing import Pool, cpu_count
import json
import time
import math
from shapely import wkb
from shapely.geometry import Point


router = APIRouter()

def get_geom(db,iszjcj,pac):
    if iszjcj=='zj':
        pac = pac[:6]
        zjcjtable = 'tp_geojson_data_qx'
    elif iszjcj=='cj':
        pac = pac[:9]
        zjcjtable = 'tp_geojson_data_zj'
    sql = f"""select name, ST_AsBinary(geometry)as geometry from {zjcjtable} where pac like '{pac}%'"""
    result = db.execute(sql)
    infos = result.fetchall()
    return  infos

@router.post("/get_info")
@router.get("/get_info")
async def get_infos(
        body = Depends(remove_xss_json),
        # zoom_level: float = Query(..., description="Zoom level for clustering"),
        # latitude_min: float = Query(..., description="Minimum latitude"),
        # latitude_max: float = Query(..., description="Maximum latitude"),
        # longitude_min: float = Query(..., description="Minimum longitude"),
        # longitude_max: float = Query(..., description="Maximum longitude"),
        # dict_value: str = Query(None),
        # option:str  = Query(None),
        db:  Session = Depends(get_db)
):
    try:
        # 根据缩放级别动态调整分组粒度
        zoom_level = float(body['zoom_level'])
        zoom_levels = {
            3: 10000,  # 全国范围
            4: 5000,
            5: 2500,
            6: 1250,
            7: 825,
            8: 412.5,
            9: 256.25,
            10: 178.125,
            11: 69.0625,
            12: 29.53125,
            13: 13.765625,
            14: 5.8828125,
            15: 2.44140625,
            16: 1.220703125,
            17: 0.6103515625,
            18: 0.30517578125
        }
        distance_threshold=zoom_levels[int(zoom_level-1)]
        # distance_threshold = 100000 / (2.2 ** zoom_level)  # 例如：每缩放一级，距离阈值减半
        dict_value= body['dict_value'].split(',')
        latitude_min = float(body['latitude_min'])
        latitude_max = float(body['latitude_max'])
        longitude_min = float(body['longitude_min'])
        longitude_max = float(body['longitude_max'])
        option = body['option'].split(',')
        print("1",time.time())
        iszjcj = ''
        pac = ''
        if 'iszjcj'in body:
            iszjcj = body['iszjcj']
            pac = body['pac']
        videos = get_videos(db,dict_value,latitude_min,latitude_max,longitude_min,longitude_max,iszjcj)
        infos = get_points(db,option,latitude_min,latitude_max,longitude_min,longitude_max,iszjcj,pac)
        # 动态分组逻辑
        if 'iszjcj' in body:
            iszjcj = body['iszjcj']
            pac = body['pac']
            groups = group_points(videos + infos, distance_threshold,db,iszjcj,pac)
        else:
            groups = group_points(videos+infos, distance_threshold)

        print("4",time.time())

        return {"code": 200,
                "msg": "操作成功",
                "data": groups}
    except Exception as e:
        # 处理异常
        traceback.print_exc()
        raise HTTPException(status_code=status.HTTP_500_INTERNAL_SERVER_ERROR, detail=str(e))
@router.post("/get_details")
@router.get("/get_details")
async def get_details(
        body = Depends(remove_xss_json),
    # center_latitude: float = Query(..., description="网格中心点的纬度"),
    # center_longitude: float = Query(..., description="网格中心点的经度"),
    # zoom_level: float = Query(..., description="缩放级别"),
    db: Session = Depends(get_db)
):
    try:
        # 计算网格大小
        zoom_level = float(body['zoom_level'])
        zoom_levels = {
            3: 10000,  # 全国范围
            4: 5000,
            5: 2500,
            6: 1250,
            7: 825,
            8: 412.5,
            9: 256.25,
            10: 178.125,
            11: 69.0625,
            12: 29.53125,
            13: 13.765625,
            14: 5.8828125,
            15: 2.44140625,
            16: 1.220703125,
            17: 0.6103515625,
            18: 0.30517578125
        }
        distance_threshold=zoom_levels[int(zoom_level-1)]
        # distance_threshold = 1000 / (1.5 ** zoom_level)  # 例如：每缩放一级，距离阈值减半
        grid_size = calculate_grid_size(distance_threshold)  # 地球半径为6371公里
        center_latitude = float(body['latitude'])
        center_longitude = float(body['longitude'])
        dict_value = body['dict_value'].split(',')
        option = body['option'].split(',')
        # 计算网格的经纬度范围
        latitude_min, latitude_max, longitude_min, longitude_max = get_grid_bounds_from_center(center_latitude, center_longitude, grid_size)
        videos = get_videos(db,dict_value,latitude_min,latitude_max,longitude_min,longitude_max)
        infos = get_points(db,option,latitude_min,latitude_max,longitude_min,longitude_max)
        return {"code": 200,
                "msg": "操作成功",
                "data": videos+infos }#{"videos":videos,"points":infos}}
    except Exception as e:
        # 处理异常
        traceback.print_exc()
        raise HTTPException(status_code=status.HTTP_500_INTERNAL_SERVER_ERROR, detail=str(e))

def calculate_grid_size(distance_threshold):
    # 假设地球半径为6371公里，将距离阈值转换为经纬度的差值
    # 这里假设纬度变化对距离的影响较小，仅根据经度计算网格大小
    earth_radius = 6371  # 地球半径，单位为公里
    grid_size = distance_threshold / earth_radius
    return grid_size

def get_grid_key(latitude, longitude, grid_size):
    # 根据经纬度和网格大小计算网格键
    return (math.floor(latitude / grid_size), math.floor(longitude / grid_size))

def get_grid_bounds_from_center(center_latitude, center_longitude, grid_size):
    half_grid_size = grid_size / 2
    min_latitude = center_latitude - half_grid_size
    max_latitude = center_latitude + half_grid_size
    min_longitude = center_longitude - half_grid_size
    max_longitude = center_longitude + half_grid_size
    return min_latitude, max_latitude, min_longitude, max_longitude

def calculate_distance(point1, point2):
    # 使用 Haversine 公式计算两点之间的距离
    from math import radians, sin, cos, sqrt, atan2
    R = 6371  # 地球半径（公里）
    lat1, lon1 = radians(point1.latitude), radians(point1.longitude)
    lat2, lon2 = radians(point2.latitude), radians(point2.longitude)
    dlat = lat2 - lat1
    dlon = lon2 - lon1
    a = sin(dlat / 2) ** 2 + cos(lat1) * cos(lat2) * sin(dlon / 2) ** 2
    c = 2 * atan2(sqrt(a), sqrt(1 - a))
    return R * c

def group_points(points, distance_threshold,db=get_db(),iszjcj='',pac=''):
    grid_size = calculate_grid_size(distance_threshold)
    grid = defaultdict(lambda:{"count":0}) #,"list":[]
    groups = []
    tmp = defaultdict(list)
    for point in points:
        grid_key = get_grid_key(float(point.latitude), float(point.longitude), grid_size)
        lovalue = str(point.latitude)+str(point.longitude)
        if lovalue not in tmp['%s-%s'%grid_key]:
            tmp['%s-%s'%grid_key].append(lovalue)
        grid['%s-%s'%grid_key]['count']+=1
        if grid['%s-%s'%grid_key]['count']>1 and len(tmp['%s-%s'%grid_key])<3:
            grid['%s-%s'%grid_key]['dataType'] = ''
            grid['%s-%s'%grid_key]['id'] = ""
            if len(tmp['%s-%s'%grid_key])<2:
                grid['%s-%s'%grid_key]['name'] = '多数据点位'
                grid['%s-%s' % grid_key]['type'] ='1'
            else:
                grid['%s-%s' % grid_key]['name'] = '聚合点位'
                grid['%s-%s' % grid_key]['type'] = '3'
            # grid['%s-%s'%grid_key]['latitude'] = float(point.latitude) #(grid_key[0] + 0.5) * grid_size
            # grid['%s-%s'%grid_key]['longitude'] = float(point.longitude) #(grid_key[1] + 0.5) * grid_size
        elif grid['%s-%s'%grid_key]['count']==1:
            if point.dataType=='video':
                grid['%s-%s' % grid_key]['id'] = point.gbIndexCode
            else:
                grid['%s-%s' % grid_key]['id'] = point.id
            grid['%s-%s'%grid_key]['dataType'] = point.dataType
            grid['%s-%s'%grid_key]['infoType'] = point.infoType
            grid['%s-%s'%grid_key]['name'] = point.name
            grid['%s-%s'%grid_key]['type'] ='2'
            grid['%s-%s'%grid_key]['latitude'] = float(point.latitude)
            grid['%s-%s'%grid_key]['longitude'] = float(point.longitude)

    groups = list(grid.values())

    if iszjcj!='':
        geom = get_geom(db, iszjcj, pac)
        for group in groups:
            for name, wkb_bytes in geom:
                lon, lat = group['longitude'],group['latitude']
                if wkb.loads(wkb_bytes).contains(Point(lon, lat)):
                    break
            else:
                groups.remove(group)
    return groups

def get_videos(db:Session,dict_value,latitude_min,latitude_max,longitude_min,longitude_max,iszjcj=''):
    que = True
    if len(dict_value)>0:
        videolist = []
        for value in dict_value:
            tag_info = get_dict_data_info(db, 'video_type', value)
            if tag_info:
                if tag_info.dict_label == '全量视频':
                    break
                else:
                    videolist += [i.video_code for i in tag_get_video_tag_list(db, value)]
        else:
            que = TPVideoInfo.gbIndexCode.in_(videolist)
    if iszjcj!='':

        pass
    # 查询分组
    query = (
        select(
            TPVideoInfo.gbIndexCode,
            TPVideoInfo.latitude,
            TPVideoInfo.longitude,
            TPVideoInfo.name,
            TPVideoInfo.status,
            literal('video').label("dataType"),
            literal('video').label("infoType")
        )
            .select_from(TPVideoInfo).where(
            and_(
                TPVideoInfo.latitude >= latitude_min,
                TPVideoInfo.latitude <= latitude_max,
                TPVideoInfo.longitude >= longitude_min,
                TPVideoInfo.longitude <= longitude_max,
                TPVideoInfo.longitude > 0,
                TPVideoInfo.latitude > 0, que
            )
        )
            .order_by(TPVideoInfo.status.asc())
    )

    result = db.execute(query)
    videos = result.fetchall()
    return videos


def get_points(db:Session,option,latitude_min,latitude_max,longitude_min,longitude_max,iszjcj='',pac=''):
    # 使用参数化查询避免 SQL 注入
    if isinstance(option, list):
        option = tuple(option)
    # if iszjcj=='zj':
    #     zd = ',T2.name as pacname,T2.pac,T2.parent_pac'
    #     pac = pac[:6]
    #     zjcjtable = f"""select * from tp_geojson_data_qx where pac like '{pac}%'"""
    #     que = f' JOIN ({zjcjtable}) T2 on ST_Intersects(T2.geometry, ST_SRID(POINT(A.longitude, A.latitude), 4326))'
    # elif iszjcj=='cj':
    #     pac = pac[:9]
    #     zjcjtable = f"""select * from tp_geojson_data_zj where pac like '{pac}%'"""
    #     que = f'LEFT JOIN ({zjcjtable}) T2 on ST_Intersects(T2.geometry, ST_SRID(POINT(A.longitude, A.latitude), 4326))'
    # else:
    que=''
    query = text(f"""
        SELECT 
            A.`name`,A.`id`,A.dataType,A.longitude,A.latitude,A.infoType 
        FROM (
            SELECT 
                *,
                ROW_NUMBER() OVER (PARTITION BY longitude, latitude, `name` 
                                   ORDER BY longitude, latitude, `name`) AS rn
            FROM 
                `point_data` 
            WHERE 
                longitude > 0 
                AND latitude BETWEEN :latitude_min AND :latitude_max 
                AND longitude BETWEEN :longitude_min AND :longitude_max 
                AND dataType IN :option 
        ) AS A  {que}
        WHERE rn = 1
    """)

    # 执行查询并传递参数
    result = db.execute(query, {
        'latitude_min': latitude_min,
        'latitude_max': latitude_max,
        'longitude_min': longitude_min,
        'longitude_max': longitude_max,
        'option': option
    })
    infos = result.fetchall()
    return infos

@router.post("/get_geojson")
async def get_geojson(
        body = Depends(remove_xss_json),
        db:  Session = Depends(get_db)
):
    try:
        # 根据缩放级别动态调整分组粒度
        latitude_min = float(body['latitude_min'])
        latitude_max = float(body['latitude_max'])
        longitude_min = float(body['longitude_min'])
        longitude_max = float(body['longitude_max'])
        if latitude_min<-90 or latitude_max>90 :
            return JSONResponse(status_code=500, content={"code": 500, "msg": "Latitude must be within [-90.000000, 90.000000]"})
        if longitude_min<-180 or longitude_max>180 :
            return JSONResponse(status_code=500, content={"code": 500, "msg": "Longitude must be within [-180.000000, 180.000000]"})
        table_name = 'tp_geojson_data_zj'
        option = body['option']
        if 'cj' == option:
            table_name = 'tp_geojson_data_cj_sq'
        sql = f"""SELECT id,
       name,
       pac,
       ST_AsGeoJSON(geometry) AS geojson,
       properties
FROM {table_name}
WHERE ST_Intersects(
        geometry,
ST_GeomFromText(
            'POLYGON(({latitude_min} {longitude_min},
                      {latitude_max} {longitude_min},
                      {latitude_max} {longitude_max},
                      {latitude_min} {longitude_max},
                      {latitude_min} {longitude_min}))',
            4326
        )
      );"""
        result = db.execute(sql)
        features = result.fetchall()

        return {"code": 200,
                "msg": "操作成功","type":"FeatureCollection",
                "features": features}
    except Exception as e:
        # 处理异常
        traceback.print_exc()
        raise HTTPException(status_code=status.HTTP_500_INTERNAL_SERVER_ERROR, detail=str(e))

@router.post("/get_geojson_new")
async def get_geojson(
        body = Depends(remove_xss_json),
        db:  Session = Depends(get_db)
):
    try:
        # 根据缩放级别动态调整分组粒度
        pac = body['area_code']
        table_name = 'tp_geojson_data_zj'
        option = body['option']
        if 'cj' == option:
            table_name = 'tp_geojson_data_cj_sq'
            pac = pac[:9]
        else:
            pac = pac[:6]
        sql = f"""SELECT 
       ST_AsGeoJSON(geometry) AS geometry,
       properties
FROM {table_name}
WHERE parent_pac = '{pac}';"""

        def gen():
            # 1. 写头
            yield '{"type":"FeatureCollection","features":['
            first = True
            # 2. 逐行流式
            for geom, prop_json in db.execute(sql):  # 迭代器，不 fetchall
                if not first:
                    yield ","
                feature = {
                    "geometry": json.loads(geom),
                    "properties": {
                        "PAC": json.loads(prop_json)['PAC'],
                        "NAME": json.loads(prop_json)['NAME']
                    }
                }
                yield json.dumps(feature, ensure_ascii=False)
                first = False
            # 3. 写尾
            yield "]}"


        return StreamingResponse(gen(), media_type="application/json",
                             headers={"Content-Disposition": "attachment; filename=data.geojson"})

    except Exception as e:
        # 处理异常
        traceback.print_exc()
        raise HTTPException(status_code=status.HTTP_500_INTERNAL_SERVER_ERROR, detail=str(e))

# @router.post("/get_map_img")
# async def get_map_img(
#     body: Dict[str, Any] = Depends(remove_xss_json),
#     db: Session = Depends(get_db),
# ):
#     """
#     输入：
#     {
#       "latitude_min": 27.0,
#       "latitude_max": 30.0,
#       "longitude_min": 118.0,
#       "longitude_max": 121.0,
#       "option": "zj"        // "zj" 或 "cj"
#     }
#     返回：PNG 图片
#     """
#     try:
#         # 1. 参数提取与合法性校验
#         lat_min = float(body["latitude_min"])
#         lat_max = float(body["latitude_max"])
#         lon_min = float(body["longitude_min"])
#         lon_max = float(body["longitude_max"])
#         option = body.get("option", "zj")
#
#         if not (-90 <= lat_min <= 90 and -90 <= lat_max <= 90):
#             raise ValueError("Latitude must be within [-90, 90]")
#         if not (-180 <= lon_min <= 180 and -180 <= lon_max <= 180):
#             raise ValueError("Longitude must be within [-180, 180]")
#
#         table_name = "tp_geojson_data_zj" if option != "cj" else "tp_geojson_data_cj_sq"
#
#         # 2. 构造 SQL
#         sql = text(
#             f"""
#             SELECT id, name, pac, ST_AsGeoJSON(geometry) AS geojson, properties
#             FROM {table_name}
#             WHERE ST_Intersects(
#                     geometry,
#                     ST_GeomFromText(
#                         'POLYGON(({lon_min} {lat_min},
#                                   {lon_max} {lat_min},
#                                   {lon_max} {lat_max},
#                                   {lon_min} {lat_max},
#                                   {lon_min} {lat_min}))',
#                         4326
#                     )
#                   );
#             """
#         )
#
#         rows = db.execute(sql).fetchall()
#         if not rows:
#             raise HTTPException(
#                 status_code=status.HTTP_404_NOT_FOUND,
#                 detail="No data within the given bbox."
#             )
#
#         # 3. 组装 GeoDataFrame
#         features = [
#             {**json.loads(r.geojson), "properties": json.loads(r.properties)}
#             for r in rows
#         ]
#         gdf = gpd.GeoDataFrame.from_features(features, crs="EPSG:4326")
#
#         # 4. 绘图
#         fig, ax = plt.subplots(figsize=(6, 6), dpi=150)
#         gdf.to_crs(epsg=3857).plot(ax=ax, alpha=0.5, edgecolor="black")
#         ctx.add_basemap(ax, source=ctx.providers.Stamen.TonerLite, crs=gdf.to_crs(epsg=3857).crs)
#         ax.set_axis_off()
#         plt.tight_layout(pad=0)
#
#         # 5. 保存成字节流
#         buf = io.BytesIO()
#         fig.savefig(buf, format="png")
#         buf.seek(0)
#         plt.close(fig)
#
#         # 6. 返回图片
#         return StreamingResponse(buf, media_type="image/png")
#
#     except Exception as e:
#         raise HTTPException(
#             status_code=status.HTTP_500_INTERNAL_SERVER_ERROR,
#             detail=str(e)
#         )

# TILE_EXTENT = 4096  # 标准 MVT 精度
# # ---------- 纯 Python 坐标转换 ----------
# def lonlat2xy(lon: float, lat: float) -> tuple[float, float]:
#     """4326 -> 3857"""
#     x = lon * 20037508.34 / 180
#     y = math.log(math.tan((90 + lat) * math.pi / 360)) / (math.pi / 180)
#     y = y * 20037508.34 / 180
#     return x, y
#
# import mercantile
# from mapbox_vector_tile import encode
# from shapely.geometry import shape
# @router.get("/tile/{option}/{z}/{x}/{y}.pbf")
# async def get_tile(
#     option: str = Path(..., regex="^(zj|cj)$"),
#     z: int = Path(..., ge=0, le=22),
#     x: int = Path(..., ge=0),
#     y: int = Path(..., ge=0),
#     db: Session = Depends(get_db),
# ):
#     """
#     根据 Slippy Map 标准 XYZ 返回 MVT 二进制。
#     前端 layer.url = '/tile/zj/{z}/{x}/{y}.pbf'
#     """
#     table = "tp_geojson_data_zj" if option == "zj" else "tp_geojson_data_cj_sq"
#
#     # 1. 计算瓦片 bbox (4326)
#     tile_bounds = mercantile.bounds(mercantile.Tile(x, y, z))
#     xmin, ymin, xmax, ymax = tile_bounds.west, tile_bounds.south, tile_bounds.east, tile_bounds.north
#     print(xmin, ymin, xmax, ymax)
#     # 2. 查相交要素
#     sql = text(
#         f"""
#         SELECT id, name, pac, properties,
#                ST_AsGeoJSON(geometry) AS geojson
#         FROM {table}
#         WHERE ST_Intersects(
#                 geometry,
#                 ST_GeomFromText(
#                     'POLYGON((
# {ymin} {xmin},
# {ymin} {xmax},
# {ymax} {xmax},
# {ymax} {xmin},
# {ymin} {xmin}))',
#                     4326
#                 )
#               );
#         """
#     )
#     rows = db.execute(sql).fetchall()
#     if not rows:
#         raise HTTPException(status_code=204)
#
#     # 3. 构造 MVT features
#     features: List[Dict[str, Any]] = []
#     bounds_3857 = mercantile.xy_bounds(mercantile.Tile(x, y, z))
#     bx, by, bw, bh = bounds_3857.left, bounds_3857.bottom, \
#                      bounds_3857.right - bounds_3857.left, \
#                      bounds_3857.top - bounds_3857.bottom
#
#     for r in rows:
#         # 直接用原始 GeoJSON，不做任何坐标变换
#         geo = json.loads(r.geojson)
#         features.append({
#             "geometry": geo,  # 必须是 4326 坐标
#             "properties": json.loads(r.properties),
#         })
#
#     # 4. 生成 MVT
#     mvt_bytes = encode([
#         {
#             "name": "layer",
#             "features": features,
#             "extent": TILE_EXTENT,
#         }
#     ])
#
#     return Response(content=mvt_bytes, media_type="application/x-protobuf")