import geopandas as gpd

# GeoParquet 的读写非常简洁
# 写入
gdf.to_parquet("data.parquet")

# 读取
gdf = gpd.read_parquet("data.parquet")

import geopandas as gpd
import pyarrow.parquet as pq

# 查看 GeoParquet 的元数据
parquet_file = pq.ParquetFile("data.parquet")
metadata = parquet_file.schema_arrow.metadata

# 解析 geo 元数据
import json
geo_metadata = json.loads(metadata[b"geo"])
print(json.dumps(geo_metadata, indent=2))

geopandas.read_parquet(
    path,               # 文件路径（本地或远程）
    columns=None,       # 要读取的列
    storage_options=None,  # 存储选项（用于远程文件）
    bbox=None,          # 边界框过滤
    **kwargs            # 传递给底层引擎的参数
)

import geopandas as gpd

# 读取本地文件
gdf = gpd.read_parquet("data.parquet")
print(f"行数: {len(gdf)}")
print(f"列数: {len(gdf.columns)}")
print(f"CRS: {gdf.crs}")
print(f"几何类型: {gdf.geom_type.unique()}")

# 只读取特定列
gdf = gpd.read_parquet(
    "large_dataset.parquet",
    columns=["name", "population", "geometry"]
)

# 对比全量读取
# 如果文件有 50 列，只读 3 列可以减少约 94% 的读取量

# 查看文件中的所有列（不读取数据）
import pyarrow.parquet as pq
schema = pq.read_schema("large_dataset.parquet")
print("可用列:")
for field in schema:
    print(f"  {field.name}: {field.type}")

# 使用边界框进行空间过滤
gdf = gpd.read_parquet(
    "china_pois.parquet",
    bbox=(115, 39, 117, 41)  # 北京周边
)
print(f"过滤后行数: {len(gdf)}")

# bbox 过滤利用 Parquet 的行组统计信息
# 可以跳过不在范围内的行组，大大提高性能

# 同时使用列选择和空间过滤
gdf = gpd.read_parquet(
    "nationwide_buildings.parquet",
    columns=["building_type", "area_sqm", "geometry"],
    bbox=(120, 30, 122, 32)  # 上海区域
)

# 这种组合使用可以极大地减少数据读取量
# 对于 TB 级别的数据集尤其有效

# Parquet 数据集可以按目录分区
# 例如：
# data/
#   province=北京/
#     part-0.parquet
#   province=上海/
#     part-0.parquet
#   province=广东/
#     part-0.parquet

# 读取整个分区数据集
gdf = gpd.read_parquet("data/")

# 使用 PyArrow 的过滤功能
import pyarrow.parquet as pq
import pyarrow.dataset as ds

dataset = ds.dataset("data/", format="parquet")
table = dataset.to_table(
    filter=ds.field("province") == "北京"
)
gdf = gpd.GeoDataFrame.from_arrow(table)

GeoDataFrame.to_parquet(
    path,                   # 输出文件路径
    index=None,             # 是否包含索引
    compression="snappy",   # 压缩算法
    schema_version=None,    # GeoParquet 规范版本
    write_covering_bbox=None,  # 是否写入覆盖边界框
    geometry_encoding=None, # 几何编码方式
    **kwargs                # 传递给底层引擎的参数
)

import geopandas as gpd
from shapely.geometry import Point

# 创建示例数据
gdf = gpd.GeoDataFrame(
    {
        "name": ["北京", "上海", "广州"],
        "population": [2171, 2487, 1868],
        "province": ["北京市", "上海市", "广东省"]
    },
    geometry=[
        Point(116.407, 39.904),
        Point(121.474, 31.230),
        Point(113.264, 23.129)
    ],
    crs="EPSG:4326"
)

# 基本写入
gdf.to_parquet("cities.parquet")

# 使用不同的压缩算法
gdf.to_parquet("data_snappy.parquet", compression="snappy")  # 默认
gdf.to_parquet("data_zstd.parquet", compression="zstd")      # 推荐
gdf.to_parquet("data_gzip.parquet", compression="gzip")      # 高压缩
gdf.to_parquet("data_none.parquet", compression=None)        # 无压缩

# 比较文件大小
import os
for comp in ["snappy", "zstd", "gzip", None]:
    filename = f"data_{comp}.parquet"
    gdf.to_parquet(filename, compression=comp)
    size = os.path.getsize(filename)
    print(f"{str(comp):10s} → {size:,} bytes")

# 不写入索引（推荐，减少文件大小）
gdf.to_parquet("data.parquet", index=False)

# 写入索引
gdf.to_parquet("data.parquet", index=True)

# 自动判断（默认行为）
gdf.to_parquet("data.parquet")  # 如果索引有意义则写入

# WKB 编码（默认，兼容性最好）
gdf.to_parquet("data_wkb.parquet", geometry_encoding="WKB")

# GeoArrow 编码（更高效，需要 GeoParquet 1.1+）
gdf.to_parquet("data_geoarrow.parquet", geometry_encoding="geoarrow")

# 写入覆盖边界框列（用于空间过滤优化）
gdf.to_parquet("data.parquet", write_covering_bbox=True)

# 这会在 Parquet 文件中额外添加 bbox 列
# 使得空间过滤查询可以利用 Parquet 的行组统计信息
# 极大地提升大数据集的空间查询性能

import geopandas as gpd

# 基本读取
gdf = gpd.read_feather("data.feather")

# 列选择
gdf = gpd.read_feather(
    "data.feather",
    columns=["name", "geometry"]
)

# 基本写入
gdf.to_feather("data.feather")

# 带压缩
gdf.to_feather("data.feather", compression="zstd")
gdf.to_feather("data.feather", compression="lz4")
gdf.to_feather("data.feather", compression="uncompressed")

# 场景1：中间结果缓存
# 在复杂的分析流程中，将中间结果保存为 Feather 以加速后续步骤
gdf_processed = complex_spatial_analysis(gdf_raw)
gdf_processed.to_feather("cache/step1_result.feather")

# 后续步骤直接读取（极快）
gdf_cached = gpd.read_feather("cache/step1_result.feather")

# 场景2：Python 进程间数据传递
# 进程 A 写入
gdf.to_feather("shared/data.feather")
# 进程 B 读取
gdf = gpd.read_feather("shared/data.feather")

import geopandas as gpd

# 使用 GeoArrow 编码写入 Parquet
gdf.to_parquet("data_geoarrow.parquet", geometry_encoding="geoarrow")

# 读取时自动识别编码方式
gdf = gpd.read_parquet("data_geoarrow.parquet")

import geopandas as gpd
import time

# 创建大数据集进行性能测试
import numpy as np
from shapely.geometry import Point

n = 1_000_000
gdf_large = gpd.GeoDataFrame(
    {"value": np.random.randn(n)},
    geometry=[Point(x, y) for x, y in
              zip(np.random.uniform(73, 135, n),
                  np.random.uniform(18, 54, n))],
    crs="EPSG:4326"
)

# WKB 编码写入
start = time.time()
gdf_large.to_parquet("test_wkb.parquet", geometry_encoding="WKB")
print(f"WKB 写入: {time.time() - start:.2f}s")

# GeoArrow 编码写入
start = time.time()
gdf_large.to_parquet("test_geoarrow.parquet", geometry_encoding="geoarrow")
print(f"GeoArrow 写入: {time.time() - start:.2f}s")

# 对比读取速度
start = time.time()
_ = gpd.read_parquet("test_wkb.parquet")
print(f"WKB 读取: {time.time() - start:.2f}s")

start = time.time()
_ = gpd.read_parquet("test_geoarrow.parquet")
print(f"GeoArrow 读取: {time.time() - start:.2f}s")

# WKB 编码将每个几何体序列化为二进制格式
from shapely import wkb
from shapely.geometry import Point

point = Point(116.4, 39.9)
wkb_bytes = wkb.dumps(point)
print(f"WKB 字节数: {len(wkb_bytes)}")
print(f"WKB 内容: {wkb_bytes.hex()}")

# 在 Parquet 中，WKB 存储为 Binary 类型的列
# 每个几何体是一个独立的二进制 blob

# 显式使用 WKB 编码（确保最大兼容性）
gdf.to_parquet("compatible.parquet", geometry_encoding="WKB")

# 显式使用 GeoArrow 编码（更高性能）
# 如果几何类型不一致，会抛出错误
try:
    gdf_mixed.to_parquet("data.parquet", geometry_encoding="geoarrow")
except ValueError as e:
    print(f"GeoArrow 编码失败: {e}")
    # 回退到 WKB
    gdf_mixed.to_parquet("data.parquet", geometry_encoding="WKB")

import geopandas as gpd
from shapely.geometry import Point

gdf = gpd.GeoDataFrame(
    {"name": ["北京", "上海"], "pop": [2171, 2487]},
    geometry=[Point(116.4, 39.9), Point(121.5, 31.2)],
    crs="EPSG:4326"
)

# 转换为 Arrow Table
table = gdf.to_arrow()
print(type(table))
# <class 'pyarrow.lib.Table'>

print(table.schema)
# name: string
# pop: int64
# geometry: binary  (WKB 编码)
# -- schema metadata --
# geo: '{"version":"1.1.0",...}'

# 使用 GeoArrow 编码
table_geoarrow = gdf.to_arrow(geometry_encoding="geoarrow")
print(table_geoarrow.schema)
# geometry 列会使用 struct 或 list 类型

import geopandas as gpd
import pyarrow as pa

# 从 Arrow Table 创建 GeoDataFrame
# （Arrow Table 必须包含 GeoParquet 兼容的 geo 元数据）
gdf = gpd.GeoDataFrame.from_arrow(table)

print(type(gdf))
# <class 'geopandas.geodataframe.GeoDataFrame'>
print(gdf.crs)
# EPSG:4326

import duckdb
import geopandas as gpd

# DuckDB 原生支持 GeoParquet
con = duckdb.connect()

# 安装并加载空间扩展
con.install_extension("spatial")
con.load_extension("spatial")

# 直接在 GeoParquet 上执行 SQL 查询
result = con.execute("""
    SELECT name, population, ST_AsText(geometry) as geom_wkt
    FROM read_parquet('cities.parquet')
    WHERE population > 10000000
""").fetchdf()

# 从 DuckDB 结果创建 GeoDataFrame
gdf = con.execute("""
    SELECT *
    FROM read_parquet('cities.parquet')
    WHERE ST_Within(
        geometry,
        ST_GeomFromText('POLYGON((115 39, 117 39, 117 41, 115 41, 115 39))')
    )
""").fetch_arrow_table()

gdf = gpd.GeoDataFrame.from_arrow(gdf)

import polars as pl
import geopandas as gpd

# Polars 可以高效读取 Parquet
df_polars = pl.read_parquet("data.parquet")

# 转换为 Arrow Table，再转为 GeoDataFrame
arrow_table = df_polars.to_arrow()
# 注意：需要确保 Arrow Table 包含正确的 geo 元数据

import geopandas as gpd
import time
import os

def benchmark_format(gdf, format_name, write_func, read_func, filepath):
    """格式性能基准测试"""
    # 写入测试
    start = time.time()
    write_func(gdf, filepath)
    write_time = time.time() - start

    # 文件大小
    if os.path.isfile(filepath):
        file_size = os.path.getsize(filepath) / (1024 * 1024)  # MB
    else:
        # 目录（如 Shapefile）
        total = sum(
            os.path.getsize(os.path.join(dirpath, f))
            for dirpath, _, filenames in os.walk(filepath)
            for f in filenames
        )
        file_size = total / (1024 * 1024)

    # 读取测试
    start = time.time()
    _ = read_func(filepath)
    read_time = time.time() - start

    print(f"{format_name:15s} | "
          f"写入: {write_time:6.2f}s | "
          f"读取: {read_time:6.2f}s | "
          f"大小: {file_size:8.2f} MB")

# 运行基准测试
# benchmark_format(gdf, "Shapefile",
#     lambda g, p: g.to_file(p),
#     lambda p: gpd.read_file(p),
#     "benchmark.shp")
# benchmark_format(gdf, "GeoJSON",
#     lambda g, p: g.to_file(p, driver="GeoJSON"),
#     lambda p: gpd.read_file(p),
#     "benchmark.geojson")
# ... 等等

# 决策辅助函数
def suggest_format(use_case):
    """根据使用场景推荐文件格式"""
    suggestions = {
        "长期存储与分析": "GeoParquet (Zstd 压缩)",
        "Web 前端展示": "GeoJSON 或 FlatGeobuf",
        "进程间数据传递": "Feather",
        "与 ArcGIS 交互": "File Geodatabase 或 Shapefile",
        "云存储与远程访问": "GeoParquet 或 FlatGeobuf",
        "与旧系统兼容": "Shapefile",
        "多图层数据": "GeoPackage",
        "大数据处理": "GeoParquet (分区)",
    }
    return suggestions.get(use_case, "GeoParquet（通用推荐）")

import geopandas as gpd

# 从 AWS S3 读取
gdf = gpd.read_parquet(
    "s3://my-bucket/data/cities.parquet",
    storage_options={
        "key": "your-access-key",
        "secret": "your-secret-key",
        "region": "us-east-1"
    }
)

# 从 Google Cloud Storage 读取
gdf = gpd.read_parquet(
    "gs://my-bucket/data/cities.parquet",
    storage_options={
        "token": "path/to/credentials.json"
    }
)

# 从 Azure Blob Storage 读取
gdf = gpd.read_parquet(
    "abfs://container@account.dfs.core.windows.net/data/cities.parquet",
    storage_options={
        "account_name": "your-account",
        "account_key": "your-key"
    }
)

# 写入到 S3
gdf.to_parquet(
    "s3://my-bucket/output/result.parquet",
    storage_options={
        "key": "your-access-key",
        "secret": "your-secret-key"
    }
)

# 从公开的 HTTP URL 读取
gdf = gpd.read_parquet(
    "https://data.example.com/open-data/buildings.parquet"
)

# 许多开放数据平台提供 GeoParquet 格式
# 例如：Overture Maps、Microsoft Building Footprints 等

import fsspec

# 使用 fsspec 的文件系统抽象
fs = fsspec.filesystem("s3", anon=True)  # 匿名访问公开数据

# 列出远程文件
files = fs.ls("s3://open-data-bucket/geo/")
print(files)

# 读取
with fs.open("s3://open-data-bucket/geo/data.parquet", "rb") as f:
    gdf = gpd.read_parquet(f)

import geopandas as gpd

# 1. 选择合适的压缩算法
# 推荐使用 Zstd — 在压缩率和速度之间取得最佳平衡
gdf.to_parquet("data.parquet", compression="zstd")

# 2. 启用覆盖边界框（用于大数据集的空间查询优化）
gdf.to_parquet("data.parquet", write_covering_bbox=True)

# 3. 优先使用 GeoArrow 编码（如果几何类型一致）
if len(gdf.geom_type.unique()) == 1:
    gdf.to_parquet("data.parquet", geometry_encoding="geoarrow")
else:
    gdf.to_parquet("data.parquet", geometry_encoding="WKB")

# 4. 不要将不必要的索引写入文件
gdf.to_parquet("data.parquet", index=False)

# 1. 只读取需要的列
gdf = gpd.read_parquet("data.parquet", columns=["name", "geometry"])

# 2. 使用空间过滤
gdf = gpd.read_parquet("data.parquet", bbox=(115, 39, 117, 41))

# 3. 对于大数据集，先查看元数据
import pyarrow.parquet as pq
pf = pq.ParquetFile("data.parquet")
print(f"行数: {pf.metadata.num_rows}")
print(f"行组数: {pf.metadata.num_row_groups}")
print(f"列数: {pf.metadata.num_columns}")

# 推荐的数据处理管道

# 1. 原始数据：可能是 Shapefile、GeoJSON 等
raw_gdf = gpd.read_file("raw_data.shp")

# 2. 数据清洗和处理
cleaned_gdf = clean_and_process(raw_gdf)

# 3. 持久化为 GeoParquet（高效存储）
cleaned_gdf.to_parquet(
    "processed/data.parquet",
    compression="zstd",
    write_covering_bbox=True
)

# 4. 后续分析直接读取 GeoParquet（极快）
gdf = gpd.read_parquet("processed/data.parquet")

# 5. 最终结果按需导出为其他格式
gdf.to_file("output/result.geojson", driver="GeoJSON")  # Web 展示
gdf.to_file("output/result.gpkg")                        # GIS 软件

def convert_to_parquet(input_path, output_path, **kwargs):
    """将任意格式转换为 GeoParquet"""
    gdf = gpd.read_file(input_path)

    # 默认选项
    defaults = {
        "compression": "zstd",
        "write_covering_bbox": True,
        "index": False
    }
    defaults.update(kwargs)

    gdf.to_parquet(output_path, **defaults)

    # 打印统计信息
    import os
    input_size = os.path.getsize(input_path) / (1024 * 1024)
    output_size = os.path.getsize(output_path) / (1024 * 1024)
    ratio = input_size / output_size if output_size > 0 else 0

    print(f"转换完成:")
    print(f"  输入: {input_path} ({input_size:.2f} MB)")
    print(f"  输出: {output_path} ({output_size:.2f} MB)")
    print(f"  压缩比: {ratio:.1f}x")
    print(f"  行数: {len(gdf)}")

# 使用示例
convert_to_parquet("old_data.shp", "new_data.parquet")