249 lines
7.5 KiB
Markdown
249 lines
7.5 KiB
Markdown
# geo_tools.core.projection 使用说明
|
||
|
||
> 坐标系查询、坐标转换、投影推荐工具，基于 [pyproj](https://pyproj4.github.io/pyproj/)。
|
||
|
||
---
|
||
|
||
## 导入方式
|
||
|
||
```python
|
||
# 推荐：从顶层包导入
|
||
from geo_tools.core.projection import (
|
||
get_crs_info,
|
||
crs_to_epsg,
|
||
transform_coordinates,
|
||
transform_point,
|
||
suggest_projected_crs,
|
||
WGS84, CGCS2000, WEB_MERCATOR, CGCS2000_UTM_50N,
|
||
)
|
||
|
||
# 或直接通过 geo_tools 导入
|
||
import geo_tools
|
||
```
|
||
|
||
---
|
||
|
||
## CRS 常量
|
||
|
||
模块内置了中国地理信息处理中最常用的 CRS 快捷常量，可直接作为参数传入所有函数：
|
||
|
||
| 常量名 | EPSG | 说明 |
|
||
|--------|------|------|
|
||
| `WGS84` | `EPSG:4326` | WGS84 地理坐标系（经纬度，最通用） |
|
||
| `CGCS2000` | `EPSG:4490` | 中国国家大地坐标系 2000（经纬度） |
|
||
| `WEB_MERCATOR` | `EPSG:3857` | Web Mercator 投影（网络地图常用，单位：米） |
|
||
| `CGCS2000_UTM_50N` | `EPSG:4508` | CGCS2000 / 3° 高斯-克吕格 50 带（单位：米） |
|
||
|
||
```python
|
||
from geo_tools.core.projection import WGS84, CGCS2000
|
||
|
||
# 直接用常量替代字符串
|
||
gdf = gdf.to_crs(CGCS2000)
|
||
```
|
||
|
||
---
|
||
|
||
## 函数说明
|
||
|
||
### `get_crs_info(crs_input)` — 查询 CRS 信息
|
||
|
||
返回坐标系的详细描述字典，方便快速了解一个未知 EPSG 的含义。
|
||
|
||
**参数**
|
||
- `crs_input`：EPSG 代码字符串（如 `"EPSG:4523"`）、整数（如 `4523`）或 proj 字符串。
|
||
|
||
**返回值**（`dict`）
|
||
|
||
| 键 | 含义 |
|
||
|----|------|
|
||
| `name` | 坐标系名称 |
|
||
| `epsg` | EPSG 整数编号（无法识别时为 `None`） |
|
||
| `unit` | 坐标单位（`degree` / `metre`） |
|
||
| `is_geographic` | 是否为地理坐标系（经纬度） |
|
||
| `is_projected` | 是否为投影坐标系（平面直角） |
|
||
| `datum` | 基准面名称 |
|
||
|
||
```python
|
||
from geo_tools.core.projection import get_crs_info
|
||
|
||
# 查询读取到的 GDB 数据的 CRS 含义
|
||
info = get_crs_info("EPSG:4523")
|
||
print(info)
|
||
# {
|
||
# 'name': 'CGCS2000 / 3-degree Gauss-Kruger zone 45',
|
||
# 'epsg': 4523,
|
||
# 'unit': 'metre',
|
||
# 'is_geographic': False,
|
||
# 'is_projected': True,
|
||
# 'datum': 'China Geodetic Coordinate System 2000'
|
||
# }
|
||
|
||
# 直接传整数
|
||
info = get_crs_info(32650)
|
||
print(info["name"]) # WGS 84 / UTM zone 50N
|
||
```
|
||
|
||
---
|
||
|
||
### `crs_to_epsg(crs_input)` — 获取 EPSG 编号
|
||
|
||
将任意 CRS 描述转为整数 EPSG 编号，无法识别时返回 `None`（不抛异常）。
|
||
|
||
```python
|
||
from geo_tools.core.projection import crs_to_epsg
|
||
|
||
epsg = crs_to_epsg("EPSG:4490")
|
||
print(epsg) # 4490
|
||
|
||
epsg = crs_to_epsg("WGS 84")
|
||
print(epsg) # 4326
|
||
|
||
epsg = crs_to_epsg("invalid_crs")
|
||
print(epsg) # None
|
||
```
|
||
|
||
---
|
||
|
||
### `transform_coordinates(xs, ys, source_crs, target_crs)` — 批量坐标转换
|
||
|
||
将一组坐标点从源坐标系批量转换到目标坐标系，返回转换后的 `(xs, ys)` 列表。
|
||
|
||
**参数**
|
||
- `xs`：X 坐标序列（地理 CRS 时为**经度**）
|
||
- `ys`：Y 坐标序列（地理 CRS 时为**纬度**）
|
||
- `source_crs`：源坐标系
|
||
- `target_crs`：目标坐标系
|
||
- `always_xy`（关键字参数）：强制按 (经度/X, 纬度/Y) 顺序处理，默认 `True`，**建议不修改**
|
||
|
||
```python
|
||
from geo_tools.core.projection import transform_coordinates, WGS84, WEB_MERCATOR
|
||
|
||
# 将北京、上海、广州的 WGS84 经纬度转为 Web Mercator 米制坐标
|
||
lons = [116.4074, 121.4737, 113.2644]
|
||
lats = [39.9042, 31.2304, 23.1291]
|
||
|
||
xs, ys = transform_coordinates(lons, lats, WGS84, WEB_MERCATOR)
|
||
print(xs) # [12959618.8, 13521606.3, 12608870.0]（单位：米）
|
||
print(ys) # [4859767.2, 3649094.2, 2641877.0]
|
||
|
||
# 国家坐标系转换：CGCS2000 经纬度 → CGCS2000 3° 高斯带（50带）
|
||
from geo_tools.core.projection import CGCS2000, CGCS2000_UTM_50N
|
||
xs_proj, ys_proj = transform_coordinates(lons, lats, CGCS2000, CGCS2000_UTM_50N)
|
||
```
|
||
|
||
---
|
||
|
||
### `transform_point(x, y, source_crs, target_crs)` — 单点坐标转换
|
||
|
||
`transform_coordinates` 的单点版本，直接返回 `(x, y)` 元组。
|
||
|
||
```python
|
||
from geo_tools.core.projection import transform_point, WGS84, CGCS2000
|
||
|
||
# 单点：WGS84 → CGCS2000（两者数值非常接近，差异在毫米级）
|
||
x, y = transform_point(116.4074, 39.9042, WGS84, CGCS2000)
|
||
print(f"CGCS2000 坐标：经度={x:.6f}, 纬度={y:.6f}")
|
||
|
||
# 单点：经纬度 → 投影坐标（米）
|
||
from geo_tools.core.projection import WEB_MERCATOR
|
||
mx, my = transform_point(116.4074, 39.9042, WGS84, WEB_MERCATOR)
|
||
print(f"墨卡托坐标：X={mx:.2f}m, Y={my:.2f}m")
|
||
```
|
||
|
||
---
|
||
|
||
### `suggest_projected_crs(lon, lat)` — 自动推荐投影 CRS
|
||
|
||
根据数据中心坐标（WGS84 经纬度）自动推荐适合**面积/距离计算**的 UTM 投影带，避免在地理坐标系下计算面积出错。
|
||
|
||
**参数**
|
||
- `lon`：中心经度（WGS84）
|
||
- `lat`：中心纬度（WGS84，北半球为正）
|
||
|
||
**返回值**：EPSG 代码字符串，如 `"EPSG:32650"`
|
||
|
||
```python
|
||
from geo_tools.core.projection import suggest_projected_crs
|
||
|
||
# 云南马关县（约 104.4°E, 23.0°N）
|
||
proj_crs = suggest_projected_crs(lon=104.4, lat=23.0)
|
||
print(proj_crs) # EPSG:32648 (WGS84 UTM zone 48N)
|
||
|
||
# 北京（116.4°E, 39.9°N）
|
||
proj_crs = suggest_projected_crs(lon=116.4, lat=39.9)
|
||
print(proj_crs) # EPSG:32650 (WGS84 UTM zone 50N)
|
||
|
||
# 实际场景：读取 GDB 后用推荐的投影计算面积
|
||
import geo_tools
|
||
|
||
gdf = geo_tools.read_gdb("data.gdb", layer="图斑")
|
||
cx, cy = gdf.geometry.unary_union.centroid.x, gdf.geometry.unary_union.centroid.y
|
||
|
||
# 如果数据是投影坐标系（单位：米），先转到地理坐标系再推荐
|
||
if gdf.crs.is_projected:
|
||
cx, cy = geo_tools.transform_point(cx, cy, gdf.crs, "EPSG:4326")
|
||
|
||
proj_crs = suggest_projected_crs(cx, cy)
|
||
gdf_proj = geo_tools.reproject(gdf, proj_crs) # 重投影
|
||
gdf_proj = geo_tools.add_area_column(gdf_proj) # 计算面积（单位：m²）
|
||
```
|
||
|
||
---
|
||
|
||
## 常见场景示例
|
||
|
||
### 场景一：不认识数据的 CRS，先查一下
|
||
|
||
```python
|
||
import geo_tools
|
||
|
||
gdf = geo_tools.read_gdb("临时数据库.gdb", layer="马关综合后图斑")
|
||
# 读取完成：CRS=EPSG:4523
|
||
|
||
info = geo_tools.get_crs_info(gdf.crs)
|
||
print(info["name"]) # CGCS2000 / 3-degree Gauss-Kruger zone 45
|
||
print(info["unit"]) # metre（投影坐标系，单位是米）
|
||
print(info["is_projected"]) # True
|
||
```
|
||
|
||
### 场景二：统一坐标系后叠置分析
|
||
|
||
```python
|
||
import geo_tools
|
||
from geo_tools.core.projection import CGCS2000
|
||
|
||
layer_a = geo_tools.read_gdb("a.gdb", layer="林地") # EPSG:4523
|
||
layer_b = geo_tools.read_vector("b.geojson") # EPSG:4326
|
||
|
||
# 统一到 CGCS2000 地理坐标系后再做叠置
|
||
layer_a = geo_tools.reproject(layer_a, CGCS2000)
|
||
layer_b = geo_tools.reproject(layer_b, CGCS2000)
|
||
|
||
result = geo_tools.overlay(layer_a, layer_b, how="intersection")
|
||
```
|
||
|
||
### 场景三：在地理坐标系数据上正确计算面积
|
||
|
||
```python
|
||
import geo_tools
|
||
from geo_tools.core.projection import suggest_projected_crs
|
||
|
||
gdf = geo_tools.read_vector("data.geojson") # EPSG:4326，单位是度
|
||
|
||
# 自动推荐合适的投影
|
||
proj = suggest_projected_crs(lon=105.0, lat=25.0) # 云贵地区
|
||
|
||
gdf = geo_tools.add_area_column(gdf, projected_crs=proj)
|
||
print(gdf["area_m2"].describe())
|
||
```
|
||
|
||
---
|
||
|
||
## 注意事项
|
||
|
||
> [!WARNING]
|
||
> 在**地理坐标系**（EPSG:4326 / 4490）下直接调用 `geometry.area` 得到的是"平方度"，**不是平方米**，面积计算会严重失真。始终用 `add_area_column()` 或先 `reproject()` 到投影坐标系后再计算。
|
||
|
||
> [!NOTE]
|
||
> `WGS84`（EPSG:4326）与 `CGCS2000`（EPSG:4490）的坐标数值差异极小（通常 < 1 米），在普通精度的分析中可视为等价，但正式国家项目中必须使用 CGCS2000。
|