使用GDAL进行重投影,配合Rasterio和Fiona#

最简单的重新投影方法是 GeoDataFrame.to_crs(). 它使用 pyproj 作为引擎,并转换几何体中的点。

这些示例演示了如何使用Fiona或rasterio作为引擎重新投影您的数据。Fiona和rasterio由GDAL提供支持,并且使用考虑几何形状的算法,而不仅仅是几何形状包含的点。这对于反经线裁剪特别有用。然而,这也意味着转换的速度不会那么快。

Fiona 示例#

from functools import partial

import fiona
import geopandas
from fiona.transform import transform_geom
from packaging import version
from pyproj import CRS
from pyproj.enums import WktVersion
from shapely.geometry import mapping, shape


# set up Fiona transformer
def crs_to_fiona(proj_crs):
    proj_crs = CRS.from_user_input(proj_crs)
    if version.parse(fiona.__gdal_version__) < version.parse("3.0.0"):
        fio_crs = proj_crs.to_wkt(WktVersion.WKT1_GDAL)
    else:
        # GDAL 3+ can use WKT2
        fio_crs = proj_crs.to_wkt()
    return fio_crs

def base_transformer(geom, src_crs, dst_crs):
    return shape(
        transform_geom(
            src_crs=crs_to_fiona(src_crs),
            dst_crs=crs_to_fiona(dst_crs),
            geom=mapping(geom),
            antimeridian_cutting=True,
        )
    )

# load natural earth land data
world = geopandas.read_file("https://naciscdn.org/naturalearth/110m/physical/ne_110m_land.zip")

destination_crs = "EPSG:3395"
forward_transformer = partial(base_transformer, src_crs=world.crs, dst_crs=destination_crs)

# Reproject to Mercator (after dropping Antartica)
world = world.drop(7)
with fiona.Env(OGR_ENABLE_PARTIAL_REPROJECTION="YES"):
    mercator_world = world.set_geometry(world.geometry.apply(forward_transformer), crs=destination_crs)

Rasterio 示例#

该示例需要 rasterio 1.2+ 和 GDAL 3+。

import geopandas
import rasterio.warp
from shapely.geometry import shape

# load example data
world = geopandas.read_file("https://naciscdn.org/naturalearth/110m/physical/ne_110m_land.zip")
# Reproject to Mercator (after dropping Antartica)
world = world.drop(7)
destination_crs = "EPSG:3395"
geometry = rasterio.warp.transform_geom(
    src_crs=world.crs,
    dst_crs=destination_crs,
    geom=world.geometry.values,
)
mercator_world = world.set_geometry(
    [shape(geom) for geom in geometry],
    crs=destination_crs,
)