import geopandas as gpd
from shapely.geometry import Point

# 创建一个地理坐标系的 GeoDataFrame
gdf = gpd.GeoDataFrame(
    {"name": ["北京", "上海"]},
    geometry=[Point(116.4, 39.9), Point(121.5, 31.2)],
    crs="EPSG:4326"  # WGS84 地理坐标系
)

print(gdf.crs)
# EPSG:4326
print(gdf.crs.is_geographic)
# True

# 将数据转换到投影坐标系
gdf_projected = gdf.to_crs("EPSG:3857")

print(gdf_projected.crs.is_projected)
# True
print(gdf_projected.crs.axis_info)
# 显示坐标轴信息，包括单位为 metre

import geopandas as gpd
from pyproj import CRS

# 读取数据
gdf = gpd.read_file("example.shp")

# 查看 CRS 类型
print(type(gdf.crs))
# <class 'pyproj.crs.crs.CRS'>

# CRS 对象包含丰富的属性
crs = gdf.crs
print(crs.name)             # 坐标系名称
print(crs.datum)            # 基准面
print(crs.ellipsoid)        # 椭球体
print(crs.is_geographic)    # 是否为地理坐标系
print(crs.is_projected)     # 是否为投影坐标系
print(crs.axis_info)        # 坐标轴信息

from pyproj import CRS

crs = CRS.from_epsg(4326)

# 基本信息
print(f"名称: {crs.name}")
# 名称: WGS 84

print(f"类型: {'地理坐标系' if crs.is_geographic else '投影坐标系'}")
# 类型: 地理坐标系

# 椭球体参数
ellipsoid = crs.ellipsoid
print(f"椭球体: {ellipsoid.name}")
print(f"长半轴: {ellipsoid.semi_major_metre} 米")
print(f"扁率倒数: {ellipsoid.inverse_flattening}")

# 基准面
datum = crs.datum
print(f"基准面: {datum.name}")

# EPSG 代码
print(f"EPSG 代码: {crs.to_epsg()}")

# 坐标轴信息
for axis in crs.axis_info:
    print(f"轴: {axis.name}, 方向: {axis.direction}, 单位: {axis.unit_name}")

from pyproj import CRS

crs1 = CRS.from_epsg(4326)
crs2 = CRS.from_epsg(3857)

# crs1 和 crs2 是不同的对象，互不影响
print(crs1 == crs2)  # False

import geopandas as gpd

# 从文件读取时，CRS 通常自动识别
gdf = gpd.read_file("china_provinces.shp")
print(gdf.crs)
# EPSG:4326

# 如果数据没有 CRS 信息
print(gdf.crs is None)
# 可能为 True，表示 CRS 未定义

from shapely.geometry import Point

# 创建没有 CRS 的 GeoDataFrame
gdf = gpd.GeoDataFrame(
    {"city": ["北京"]},
    geometry=[Point(116.4, 39.9)]
)
print(gdf.crs)
# None

# 声明坐标系为 WGS84
gdf = gdf.set_crs("EPSG:4326")
print(gdf.crs)
# EPSG:4326

# 注意：坐标值没有变化！
print(gdf.geometry.iloc[0])
# POINT (116.4 39.9)

# 基本用法：声明 CRS
gdf = gdf.set_crs("EPSG:4326")

# 如果已有 CRS，默认会报错
# gdf = gdf.set_crs("EPSG:3857")  # 会抛出 ValueError

# 使用 allow_override=True 强制覆盖已有 CRS
gdf = gdf.set_crs("EPSG:3857", allow_override=True)

# 原地修改（修改原对象，不创建新对象）
gdf.set_crs("EPSG:4326", inplace=True, allow_override=True)

# 错误示范：用 set_crs() 企图进行坐标转换
# 这只会修改 CRS 标签，但坐标值不变，导致数据位置完全错误！
gdf_wrong = gdf.set_crs("EPSG:3857", allow_override=True)
# 坐标仍然是 (116.4, 39.9)，但 GeoPandas 会认为这是 Web Mercator 的坐标

# 正确做法：用 to_crs() 进行坐标转换
gdf_correct = gdf.to_crs("EPSG:3857")
# 坐标被正确转换为 Web Mercator 坐标

import geopandas as gpd
from shapely.geometry import Point

# 创建 WGS84 坐标数据
cities = gpd.GeoDataFrame(
    {
        "city": ["北京", "上海", "广州", "成都"],
        "population": [2171, 2487, 1868, 2094]
    },
    geometry=[
        Point(116.407, 39.904),
        Point(121.474, 31.230),
        Point(113.264, 23.129),
        Point(104.066, 30.573)
    ],
    crs="EPSG:4326"
)

print("原始坐标（WGS84）:")
print(cities.geometry.values)

# 转换到 Web Mercator
cities_mercator = cities.to_crs("EPSG:3857")
print("\n转换后坐标（Web Mercator）:")
print(cities_mercator.geometry.values)

# 适用于 Google Maps、OpenStreetMap 等 Web 地图
gdf_web = gdf.to_crs(epsg=3857)

# 北京位于 UTM 第 50 带北半球
gdf_utm = gdf.to_crs(epsg=32650)

# CGCS2000 与 WGS84 非常接近，但不完全相同
gdf_wgs84 = gdf_cgcs2000.to_crs(epsg=4326)

# 使用自定义的兰伯特等角圆锥投影（适合中国区域）
proj_string = "+proj=lcc +lat_1=25 +lat_2=47 +lat_0=36 +lon_0=105 +datum=WGS84 +units=m"
gdf_lcc = gdf.to_crs(proj_string)

# 使用 EPSG 代码
gdf.to_crs(epsg=4326)

# 使用 CRS 字符串
gdf.to_crs("EPSG:4326")

# 使用 pyproj.CRS 对象
from pyproj import CRS
target_crs = CRS.from_epsg(4326)
gdf.to_crs(target_crs)

# 使用 WKT 字符串
gdf.to_crs(wkt_string)

# 使用 Proj 字符串
gdf.to_crs("+proj=utm +zone=50 +datum=WGS84 +units=m")

# 1. 必须先有 CRS 才能转换
gdf_no_crs = gpd.GeoDataFrame(geometry=[Point(0, 0)])
# gdf_no_crs.to_crs("EPSG:3857")  # 报错！因为不知道原始 CRS

# 正确做法：先声明再转换
gdf_no_crs = gdf_no_crs.set_crs("EPSG:4326")
gdf_converted = gdf_no_crs.to_crs("EPSG:3857")

# 2. 转换可能引入微小的数值误差
# 这是正常现象，因为涉及复杂的数学变换

# 3. 某些投影会导致特定区域的较大变形
# 例如 Web Mercator 在极地区域变形严重

from pyproj import CRS

crs_4326 = CRS.from_epsg(4326)
print(crs_4326.name)
# WGS 84
print(crs_4326.ellipsoid.name)
# WGS 84
print(crs_4326.datum.name)
# World Geodetic System 1984

crs_3857 = CRS.from_epsg(3857)
print(crs_3857.name)
# WGS 84 / Pseudo-Mercator

# Web Mercator 的特点：
# - 保角性好（形状正确）
# - 面积变形严重（尤其在高纬度地区）
# - 不适合用于面积计算

crs_4490 = CRS.from_epsg(4490)
print(crs_4490.name)
# China Geodetic Coordinate System 2000
print(crs_4490.datum.name)
# China 2000

# CGCS2000 与 WGS84 的椭球体参数非常接近
# WGS84 长半轴: 6378137.0 m, 扁率倒数: 298.257223563
# CGCS2000 长半轴: 6378137.0 m, 扁率倒数: 298.257222101
# 差异在毫米级别

# 根据城市位置选择对应的 UTM 分区
# 北京（经度 116.4°）位于第 50 带
gdf_beijing_utm = gdf_beijing.to_crs(epsg=32650)

# 上海（经度 121.5°）位于第 51 带
gdf_shanghai_utm = gdf_shanghai.to_crs(epsg=32651)

import geopandas as gpd
from shapely.geometry import Point

# 创建北京区域的数据
gdf = gpd.GeoDataFrame(
    {"name": ["天安门", "故宫", "颐和园"]},
    geometry=[
        Point(116.391, 39.907),
        Point(116.397, 39.917),
        Point(116.275, 39.999)
    ],
    crs="EPSG:4326"
)

# 自动估算最合适的 UTM CRS
utm_crs = gdf.estimate_utm_crs()
print(utm_crs)
# EPSG:32650（UTM zone 50N）

# 直接用估算的 CRS 进行转换
gdf_utm = gdf.to_crs(utm_crs)
print(gdf_utm.crs)
# EPSG:32650

# 基本用法
utm_crs = gdf.estimate_utm_crs()

# 指定基准面（默认为 WGS 84）
utm_crs_nad83 = gdf.estimate_utm_crs(datum_name="NAD83")

# 该方法通过计算数据的中心经度来确定 UTM 分区
# 公式：zone = floor((中心经度 + 180) / 6) + 1

# 场景1：计算面积前自动选择投影
def calculate_area_in_sqm(gdf):
    """计算面积（平方米），自动选择合适的 UTM 投影"""
    if gdf.crs.is_geographic:
        utm_crs = gdf.estimate_utm_crs()
        gdf_projected = gdf.to_crs(utm_crs)
    else:
        gdf_projected = gdf
    return gdf_projected.area

# 场景2：为不同地区的数据分别估算最佳投影
regions = {
    "华北": gdf_north,
    "华东": gdf_east,
    "西南": gdf_southwest
}
for name, region_gdf in regions.items():
    utm = region_gdf.estimate_utm_crs()
    print(f"{name}区域最佳 UTM: {utm}")

# 字符串形式
gdf.set_crs("EPSG:4326")
gdf.to_crs("EPSG:4326")

# 整数形式（通过 epsg 参数）
gdf.set_crs(epsg=4326)
gdf.to_crs(epsg=4326)

# 通过 pyproj.CRS
from pyproj import CRS
crs = CRS.from_epsg(4326)
gdf.set_crs(crs)

from pyproj import CRS

crs = CRS.from_epsg(4326)
wkt = crs.to_wkt()
print(wkt)

# 使用 WKT 字符串创建 CRS
crs_from_wkt = CRS.from_wkt(wkt)

# 也可以获取旧版 WKT1 格式
wkt1 = crs.to_wkt(version="WKT1_GDAL")

from pyproj import CRS

# 从 Proj 字符串创建 CRS
crs = CRS.from_proj4("+proj=longlat +datum=WGS84 +no_defs")

# 获取 Proj 字符串
proj_str = crs.to_proj4()
print(proj_str)
# +proj=longlat +datum=WGS84 +no_defs +type=crs

# 自定义投影
custom_proj = CRS.from_proj4(
    "+proj=tmerc +lat_0=0 +lon_0=117 +k=1 +x_0=500000 +y_0=0 "
    "+ellps=GRS80 +units=m +no_defs"
)

from pyproj import CRS

# EPSG authority
crs = CRS.from_authority("EPSG", 4326)

# OGC authority
crs_ogc = CRS.from_authority("OGC", "CRS84")

# 获取 authority 信息
auth = crs.to_authority()
print(auth)
# ('EPSG', '4326')

from pyproj import CRS

# EPSG:4326 的轴序是 纬度(Lat), 经度(Lon)
crs_4326 = CRS.from_epsg(4326)
for axis in crs_4326.axis_info:
    print(f"{axis.name}: {axis.direction}")
# Geodetic latitude: north
# Geodetic longitude: east

# OGC:CRS84 的轴序是 经度(Lon), 纬度(Lat)
crs84 = CRS.from_authority("OGC", "CRS84")
for axis in crs84.axis_info:
    print(f"{axis.name}: {axis.direction}")
# Geodetic longitude: east
# Geodetic latitude: north

# 在实际应用中，GeoPandas 始终使用 (x, y) 即 (经度, 纬度) 的顺序
# 所以在 GeoPandas 中两者的行为基本一致

from pyproj import CRS

crs1 = CRS.from_epsg(4326)
crs2 = CRS.from_proj4("+proj=longlat +datum=WGS84 +no_defs")
crs3 = CRS.from_epsg(3857)

# 直接比较
print(crs1 == crs2)  # True（虽然创建方式不同，但本质相同）
print(crs1 == crs3)  # False

# 精确比较与宽松比较
print(crs1.equals(crs2))  # True

import geopandas as gpd
from shapely.geometry import Point, Polygon

# 创建两个不同 CRS 的 GeoDataFrame
gdf1 = gpd.GeoDataFrame(
    geometry=[Point(116.4, 39.9)],
    crs="EPSG:4326"
)
gdf2 = gpd.GeoDataFrame(
    geometry=[Polygon([(0, 0), (1, 0), (1, 1), (0, 1)])],
    crs="EPSG:3857"
)

# 尝试执行空间连接 → 会抛出警告或错误
# gpd.sjoin(gdf1, gdf2)
# UserWarning: CRS mismatch between the CRS of left geometries
# and the CRS of right geometries.

# 正确做法：统一 CRS 后再操作
gdf2_4326 = gdf2.to_crs("EPSG:4326")
result = gpd.sjoin(gdf1, gdf2_4326)

# 当 CRS 为 None 时，某些操作会发出警告
gdf_no_crs = gpd.GeoDataFrame(geometry=[Point(0, 0)])
print(gdf_no_crs.crs)
# None

# 与有 CRS 的数据进行操作时会警告
# 建议始终为数据设置正确的 CRS

def validate_crs(gdf, expected_epsg=None):
    """验证 GeoDataFrame 的 CRS"""
    if gdf.crs is None:
        print("⚠️ 警告：CRS 未定义")
        return False

    print(f"✓ CRS: {gdf.crs.name}")
    print(f"  EPSG: {gdf.crs.to_epsg()}")
    print(f"  类型: {'地理坐标系' if gdf.crs.is_geographic else '投影坐标系'}")

    if expected_epsg and gdf.crs.to_epsg() != expected_epsg:
        print(f"⚠️ 期望 EPSG:{expected_epsg}，实际为 EPSG:{gdf.crs.to_epsg()}")
        return False

    return True

# 使用示例
validate_crs(gdf, expected_epsg=4326)

from pyproj import CRS

# CGCS2000 地理坐标系
cgcs2000 = CRS.from_epsg(4490)
print(cgcs2000.name)
# China Geodetic Coordinate System 2000

# CGCS2000 投影坐标系（高斯-克吕格 3° 带）
# 以中央经线 117° 为例
cgcs2000_proj = CRS.from_epsg(4548)  # CGCS2000 / 3-degree Gauss-Kruger zone 39
print(cgcs2000_proj.name)

# 北京 54 地理坐标系
# 注意：pyproj 中的 EPSG 代码可能不完整
# 可能需要使用 WKT 或 Proj 字符串定义

# 使用 Proj 字符串
bj54_proj = "+proj=longlat +ellps=krass +no_defs"

from pyproj import CRS

# 西安 80 地理坐标系
xian80 = CRS.from_epsg(4610)
print(xian80.name)
# Xian 1980

from pyproj import CRS

# CGCS2000 高斯-克吕格 6° 带第 20 带（中央经线 117°）
gk_6degree = CRS.from_epsg(4498)

# CGCS2000 高斯-克吕格 3° 带第 39 带（中央经线 117°）
gk_3degree = CRS.from_epsg(4549)

# 如果 EPSG 代码不够用，可以自定义
custom_gk = CRS.from_proj4(
    "+proj=tmerc +lat_0=0 +lon_0=117 +k=1.0 "
    "+x_0=39500000 +y_0=0 +ellps=GRS80 +units=m +no_defs"
)

import geopandas as gpd
from shapely.geometry import Point

# 创建 CGCS2000 数据
gdf_cgcs = gpd.GeoDataFrame(
    {"name": ["测试点"]},
    geometry=[Point(116.4, 39.9)],
    crs="EPSG:4490"
)

# CGCS2000 → WGS84
# 由于两个椭球体参数非常接近，差异在毫米级别
gdf_wgs84 = gdf_cgcs.to_crs("EPSG:4326")

# CGCS2000 地理坐标 → CGCS2000 投影坐标
gdf_proj = gdf_cgcs.to_crs("EPSG:4549")  # 3° 带 39 带

# 西安 80 → CGCS2000
# 需要注意：这种转换可能需要地方性的转换参数
gdf_xian80 = gpd.GeoDataFrame(
    geometry=[Point(116.4, 39.9)],
    crs="EPSG:4610"
)
gdf_to_cgcs = gdf_xian80.to_crs("EPSG:4490")

import geopandas as gpd

# ❌ 错误：直接在地理坐标系下计算面积
gdf = gpd.read_file("polygons.shp")  # CRS: EPSG:4326
area_wrong = gdf.area  # 结果单位是"平方度"，没有实际意义！

# ✅ 正确：先转换到投影坐标系再计算
utm_crs = gdf.estimate_utm_crs()
gdf_utm = gdf.to_crs(utm_crs)
area_correct = gdf_utm.area  # 结果单位是平方米
print(f"面积: {area_correct.iloc[0]:.2f} 平方米")

# 常见原因：两个数据集的 CRS 不一致
gdf1 = gpd.read_file("data1.shp")  # EPSG:4326
gdf2 = gpd.read_file("data2.shp")  # EPSG:4490

# 检查 CRS
print(f"数据集1: {gdf1.crs}")
print(f"数据集2: {gdf2.crs}")

# 统一 CRS
gdf2 = gdf2.to_crs(gdf1.crs)

# 可能原因1：使用了错误的 CRS 标签
# 例如数据实际是 CGCS2000，但标记为 WGS84
# 检查方法：查看坐标值是否在合理范围内

# 可能原因2：投影坐标被误标为地理坐标
# 例如坐标值为 (39500116, 4430000)，这明显是投影坐标而非经纬度

# 可能原因3：经纬度顺序颠倒
# 标准经纬度：经度 73°-135°，纬度 18°-54°（中国范围）

# 读取数据后第一步就检查 CRS
gdf = gpd.read_file("data.shp")
print(f"CRS: {gdf.crs}")
print(f"EPSG: {gdf.crs.to_epsg() if gdf.crs else 'None'}")

# 存储和交换数据时使用 WGS84 或 CGCS2000
gdf.to_crs("EPSG:4326").to_file("output.geojson", driver="GeoJSON")

# 进行空间分析时转换到投影坐标系
gdf_analysis = gdf.to_crs(gdf.estimate_utm_crs())

# 定义统一的 CRS
TARGET_CRS = "EPSG:4326"

# 所有数据集统一转换
datasets = [gdf1, gdf2, gdf3]
datasets_unified = [ds.to_crs(TARGET_CRS) for ds in datasets]

# ✅ 推荐
gdf.to_crs("EPSG:4326")

# ⚠️ 不推荐（可能丢失信息）
gdf.to_crs("+proj=longlat +datum=WGS84 +no_defs")

from pyproj import CRS

crs = CRS.from_epsg(32650)
area_of_use = crs.area_of_use
print(f"适用范围: {area_of_use.name}")
print(f"经度: {area_of_use.west}° - {area_of_use.east}°")
print(f"纬度: {area_of_use.south}° - {area_of_use.north}°")