我用Python把一张省级TIF地图切成了1000张小图，发现了一些有意思的事

故事是这样的。

前两天有个朋友找我，说他手头有一张省级的高程TIF数据，几十G，想按行政区划切出来，每个区县单独一张图。我说这还不简单，用ArcGIS或者QGIS的裁剪工具不就行了？

他说，不行，要切100多个区县，手动点100多次？

我：。。。

行吧，那就写个脚本呗。

1. 先说说我们要干啥

这事儿的核心逻辑其实挺简单的，就三步：

1. 读取GeoJSON，拿到每个区县的边界坐标
2. 用这个边界去裁剪TIF影像
3. 保存成单独的图片

但真写起来，坑还是不少的。我跟你说。

2. 读取GeoJSON，遍历每个行政区划

GeoJSON这玩意儿，标准格式，长这样：

with open('./Data/quhua/xxx.geojson', 'r', encoding='utf-8') as f:
    geojson_data = json.load(f)

# 遍历 FeatureCollection 的 features
for feature in geojson_data.get("features", []):
    properties = feature.get("properties", {})
    geometry = feature.get("geometry", {})

# 获取属性信息
    adcode = properties.get("adcode")
    name = properties.get("name")
    level = properties.get("level")

# 获取边界信息
    boundaries = geometry.get("coordinates")
    geom = shape(geometry)
    bounds = geom.bounds

print(f"行政区划: {name}, 编码: {adcode}, 等级: {level}")

这里有个需要注意的点，GeoJSON里的geometry类型可能是Polygon（单多边形）也可能是MultiPolygon（多个多边形组成）。得分别处理：

if geometry.get("type") == "Polygon":
    polygon = Polygon(boundaries[0])  # 外环是第一个元素
elif geometry.get("type") == "MultiPolygon":
    polygon = MultiPolygon([Polygon(coords[0]) for coords in boundaries])

我一开始没注意这个，结果跑到某个多岛屿的区县直接报错，给我一下子整不会了。

3. 用rasterio裁剪TIF

这块是核心代码，用rasterio的mask函数：

tif_Dem_path = "/path/to/your/large.tif"

with rasterio.open(tif_Dem_path) as src:
# 使用 GeoJSON 边界掩膜裁剪
    out_image, out_transform = mask(src, [geometry], crop=True)

# 保存路径处理
prefixCode = str(adcode)[:4]
newavepath = savepath + "\\" + prefixCode
saveDemPath = newavepath + "\\dem\\"

mask函数会自动根据你给的geometry去裁剪，crop=True表示只保留有数据的部分。out_transform是新的仿射变换参数，后面画图的时候要用到。

4. 可视化并保存

光裁剪出来还不够，得画成图。我用的matplotlib + cartopy：

def drawDemByPolygon(tifname, data, bounds, transform, polygon, savePath):
    tifname = tifname + "dem"
    data = data[0]  # 取出第一波段

# 自定义颜色映射
    rainbow_cmap = plt.cm.rainbow
    colors = rainbow_cmap(np.linspace(0, 1, 256))
    colors[0] = [0, 0, 0, 0]  # 把0值设成透明
    my_cmap = ListedColormap(colors)

# 创建图形，使用PlateCarree投影
    fig, ax = plt.subplots(subplot_kw={'projection': ccrs.PlateCarree()})

# 绘制数据
    c = ax.imshow(data, 
                  extent=[transform[2], transform[2] + transform[0] * data.shape[1],
                          transform[5] + transform[4] * data.shape[0], transform[5]], 
                  cmap=my_cmap)

# 叠加边界线
if polygon:
        feature = ShapelyFeature(polygon, ccrs.PlateCarree(), 
                                edgecolor='black', facecolor='none', linewidth=1)
        ax.add_feature(feature)

# 添加colorbar
    cbar = plt.colorbar(c, ax=ax, orientation='vertical', pad=0.02)
    cbar.set_label("高程图例")

    plt.axis('off')

# 保存
    plt.savefig(f"{savePath}{tifname}.png",
                dpi=500,
                facecolor='white',
                bbox_inches='tight',
                pad_inches=0.0)

    plt.close()
    gc.collect()  # 手动回收内存，大TIF容易爆

这里我踩了个坑，一开始没加gc.collect()，跑到一半内存炸了。后来每次画完强制回收一下，稳得很。

5. 一个小技巧：扩展边界

有时候你只想看某个区县，但直接裁剪出来的图边界太贴边了，不好看。可以加个扩展边界的函数：

def expand_boundary(geometry, bounds, expansion_ratio=0.1):
    lonmin, latmin, lonmax, latmax = bounds

# 根据比例扩展
    lon_range = lonmax - lonmin
    lat_range = latmax - latmin

    lonmin_expanded = lonmin - lon_range * expansion_ratio
    lonmax_expanded = lonmax + lon_range * expansion_ratio
    latmin_expanded = latmin - lat_range * expansion_ratio
    latmax_expanded = latmax + lat_range * expansion_ratio

# 创建扩展后的矩形
    expanded_geometry = Polygon([
        (lonmin_expanded, latmin_expanded),
        (lonmax_expanded, latmin_expanded),
        (lonmax_expanded, latmax_expanded),
        (lonmin_expanded, latmax_expanded)
    ])

return expanded_geometry

这样切出来的图会带一点「呼吸感」，看着舒服多了。

6. 完整流程串起来

# 1. 读取GeoJSON
with open(geojson_path, 'r') as f:
    geojson_data = json.load(f)

# 2. 遍历每个区县
for feature in geojson_data.get("features", []):
    properties = feature.get("properties", {})
    geometry = feature.get("geometry", {})

    name = properties.get("name")
    adcode = properties.get("adcode")

# 3. 可选：扩展边界
    geom = shape(geometry)
    bounds = geom.bounds
    expanded_geometry = expand_boundary(geometry, bounds, expansion_ratio=0.1)

# 4. 裁剪TIF
with rasterio.open(tif_path) as src:
        out_image, out_transform = mask(src, [expanded_geometry], crop=True)

# 5. 画图保存
    drawDemByPolygon(name, out_image, bounds, out_transform, geom, save_path)

print(f"{name} 处理完成！")

跑一遍，100多个区县自动切完，喝杯茶的功夫。

7. 一些碎碎念

说实话，这种活儿用ArcGIS也能做，但写脚本的好处是啥呢？

一是可复现，这次跑完下次换个省接着跑。二是可定制，想加什么功能自己改，比如我后来还加了批量加水印、自动压缩、生成缩略图这些。三是。。。爽啊，看着控制台一行行输出「某某区县处理完成」，那种感觉太爽了。

而且我发现，很多看似要用专业软件才能干的活儿，其实用Python几行代码就能搞定。GIS这行尤其明显，ArcGIS、QGIS能干的事，gdal、rasterio、geopandas基本都能干，还更灵活。

当然，不是说专业软件不好用，而是有时候你只是想快速切几张图，打开ArcGIS等它启动的时间，脚本都跑完了。

效率这事儿，有时候就是这么一点点抠出来的。

最后

代码我整理了一下，需要的可以自取。主要依赖就这几个：

rasterio
matplotlib
cartopy
shapely
numpy

安装完直接跑，改一下路径就行。

以上，既然看到这里了，如果觉得不错，随手点个赞、在看、转发三连吧，如果想第一时间收到推送，也可以给我个星标⭐～

谢谢你看我的文章，我们，下次再见。

/ 作者：史蒂夫AI
/ 投稿或爆料，请联系公众号：史蒂夫AI

结语 可批量处理,个例展示 1.无扩展边界

2.扩展边界