""" geo_tools.core.projection ~~~~~~~~~~~~~~~~~~~~~~~~~ 坐标系与投影转换工具,基于 pyproj。 """ from __future__ import annotations from typing import Sequence from pyproj import CRS, Transformer import geopandas as gpd from app.utils.logger import get_logger logger = get_logger(__name__) def get_crs_info(crs_input: str | int) -> dict[str, str | int | None]: """获取 CRS 的基本信息。 Parameters ---------- crs_input: EPSG 代码(整数或 ``"EPSG:4326"`` 字符串)或 proj 字符串。 Returns ------- dict 包含 ``name``、``epsg``、``unit``、``is_geographic``、``is_projected``、``datum``。 """ crs = CRS.from_user_input(crs_input) return { "name": crs.name, "epsg": crs.to_epsg(), "unit": str(crs.axis_info[0].unit_name) if crs.axis_info else None, "is_geographic": crs.is_geographic, "is_projected": crs.is_projected, "datum": crs.datum.name if crs.datum else None, } def crs_to_epsg(crs_input: str | int) -> int | None: """尝试将 CRS 转为 EPSG 整数编号,无法识别时返回 None。""" try: return CRS.from_user_input(crs_input).to_epsg() except Exception: return None def transform_coordinates( xs: Sequence[float], ys: Sequence[float], source_crs: str | int, target_crs: str | int, *, always_xy: bool = True, ) -> tuple[list[float], list[float]]: """批量转换坐标点。 Parameters ---------- xs: X 坐标序列(地理 CRS 时为经度)。 ys: Y 坐标序列(地理 CRS 时为纬度)。 source_crs: 源 CRS。 target_crs: 目标 CRS。 always_xy: 强制以 (X, Y) 顺序输入输出(推荐保持 True)。 Returns ------- (list[float], list[float]) 转换后的 (xs, ys)。 """ transformer = Transformer.from_crs(source_crs, target_crs, always_xy=always_xy) result_xs, result_ys = transformer.transform(list(xs), list(ys)) return list(result_xs), list(result_ys) def transform_point( x: float, y: float, source_crs: str | int, target_crs: str | int, *, always_xy: bool = True, ) -> tuple[float, float]: """转换单个坐标点。""" xs, ys = transform_coordinates([x], [y], source_crs, target_crs, always_xy=always_xy) return xs[0], ys[0] def suggest_projected_crs(lon: float, lat: float, use_3degree: bool = True) -> str: """根据经纬度范围自动推荐适合面积/距离计算的投影 CRS(CGCS2000 高斯-克吕格 带号)。 Parameters ---------- lon: 中心经度(CGCS2000)。 lat: 中心纬度(CGCS2000)。 use_3degree: True 表示3度分带,False 表示6度分带。 Returns ------- str EPSG 代码字符串,如 ``"EPSG:32650"``(CGCS2000 高斯-克吕格 带号)。 """ if use_3degree: # 3度分带计算:中央经线 = 3° * n central_meridian = round(lon / 3) * 3 zone_number = int(central_meridian / 3) # CGCS2000 3度带投影定义 # 从第25带到45带(75°E-135°E) if 75 <= central_meridian <= 135: epsg = 4513 + zone_number - 25 else: # 默认使用36带(108°E) epsg = 4524 logger.warning("经度范围超出3度带范围,默认使用36带(108°E)") else: # 6度分带计算:中央经线 = 6° * n - 3° central_meridian = round((lon + 3) / 6) * 6 - 3 zone_number = int((central_meridian + 3) / 6) # CGCS2000 6度带投影定义 # 从第13带到23带(75°E-135°E) if 75 <= central_meridian <= 135: epsg = 4491 + zone_number - 13 else: # 默认使用18带(105°E) epsg = 4496 logger.warning("经度范围超出6度带范围,默认使用18带(105°E)") logger.debug("建议投影 CRS:EPSG:%d(lon=%.2f, lat=%.2f)", epsg, lon, lat) return f"EPSG:{epsg}" def reproject_gdf( gdf: gpd.GeoDataFrame, target_crs: str | int | None = None, *, auto_crs: bool = False, verbose: bool = True, ) -> gpd.GeoDataFrame: """将 GeoDataFrame(要素类)重投影到目标坐标系。 Parameters ---------- gdf: 输入 GeoDataFrame,必须已定义 CRS。 target_crs: 目标 CRS,如 ``"EPSG:4326"``、``"EPSG:4490"`` 或整数 ``4523``。 与 ``auto_crs=True`` 二选一。 auto_crs: 为 ``True`` 时忽略 ``target_crs``,根据数据中心点自动推荐 CGCS2000 高斯-克吕格 带号(默认使用 3度分带)。 verbose: 为 ``True`` 时在日志中打印投影前后的 CRS 信息。 Returns ------- gpd.GeoDataFrame 重投影后的新 GeoDataFrame(原始对象不变)。 Raises ------ ValueError ``gdf`` 未定义 CRS,或 ``target_crs`` 与 ``auto_crs`` 均未指定。 Examples -------- >>> # 指定目标 CRS >>> gdf_proj = reproject_gdf(gdf, "EPSG:4490") >>> # 自动推荐 CGCS2000 高斯-克吕格 带号 默认使用 3度分带 >>> gdf_utm = reproject_gdf(gdf, auto_crs=True) >>> # 配合 GDB 读取 >>> gdf = geo_tools.readers.read_vector("data.gdb/图斑") >>> gdf_proj = reproject_gdf(gdf, "EPSG:4326") """ if gdf.crs is None: raise ValueError("GeoDataFrame 未定义 CRS,请先调用 set_crs() 设置坐标系。") if auto_crs: # 先统一到地理坐标系,再取中心点推荐 CGCS2000 高斯-克吕格 带号 if gdf.crs.is_projected: center = gdf.to_crs("EPSG:4490").geometry.union_all().centroid else: center = gdf.geometry.union_all().centroid target_crs = suggest_projected_crs(center.x, center.y) logger.info("auto_crs:自动推荐投影 CRS = %s", target_crs) if target_crs is None: raise ValueError("请指定 target_crs,或设置 auto_crs=True 自动推荐投影。") src_crs_str = gdf.crs.to_string() result = gdf.to_crs(target_crs) if verbose: tgt_info = get_crs_info(target_crs) logger.info( "要素类重投影完成:%s → %s(%s,单位:%s,要素数:%d)", src_crs_str, tgt_info.get("epsg") or target_crs, tgt_info.get("name"), tgt_info.get("unit"), len(result), ) return result