企业官网建设流程全解析

简介：这份中国自然保护区空间分布数据以shp矢量格式整理，面向从事GIS制图、生态地理研究与空间分析的学生、科研人员及规划从业者，可用于自然保护区分布格局展示、图层叠加分析与专题地图制作。资源包共16个文件，约7.95MB，涵盖国家级与地方级两类保护区数据，主要文件类型包括shp几何数据、dbf属性表、shx索引、prj投影信息，以及cpg编码、sbn与sbx空间索引、xml元数据等配套文件，各类型协同保证数据在ArcGIS等平台中正常读取、渲染与查询。目前已有948人学习下载，说明其在同类数据中具备一定参考价值。读者可据此快速获取现成的保护区边界图层，省去繁琐的数字化采集环节，直接用于全国或区域尺度的空间分布制图、缓冲区分析与生态保护研究，也可作为教学演示与课程作业的基础数据，帮助把精力集中在分析思路与制图表达上。

1. 中国自然保护区空间分布数据shp：从一份边界文件到可复现的空间分析

手里拿到一份「中国自然保护区空间分布数据shp」，很多人第一反应是打开 ArcGIS 或 QGIS 看一眼边界长什么样，然后就没有然后了。真正卡住人的从来不是「有没有数据」，而是这份 shp 到底能不能直接进分析流程：坐标系是地理坐标还是投影坐标、属性表里有没有保护区名称和级别、边界是面还是多部件面、跨省保护区怎么处理。我见过太多人把 shp 拖进软件，出图挺好看，一到做叠加分析、面积统计、渔网分割就各种报错，最后发现是坐标系没统一或者几何有自相交。这篇笔记就围绕这份数据，把从加载、检查、坐标转换到实际空间统计的完整链路讲清楚，顺带把 shp 转 txt、shp 导出 wkt、渔网分割这些高频操作串进去。适合做生态评估、国土空间规划、保护区选址、遥感变化检测的从业者，也适合刚接触矢量数据、想拿真实数据练手的新手。读完你应该能独立判断一份保护区 shp 能不能用、该怎么用、哪里会翻车。

2. 先搞懂这份 shp 里到底装了什么：字段、几何与坐标系的三角关系

2.1 保护区 shp 的典型字段结构和几何类型

中国自然保护区空间分布数据通常以面要素（Polygon）为主，少数以点要素记录中心点。面要素的属性表一般包含这几类字段：保护区名称（NAME 或 PA_NAME）、级别（国家级/省级/市县级，字段名常见 LEVEL 或 GRADE）、类型（森林生态、野生动物、湿地、荒漠、海洋海岸等，字段名 TYPE 或 CATEGORY）、所属省份（PROVINCE）、批复面积（AREA 或 SHAPE_AREA）、建立年份。不同来源的数据字段名差异很大，有的用拼音缩写，有的用中文，有的干脆只有 OBJECTID 和几何。拿到数据第一件事不是画图，是打开属性表看字段，确认有没有能唯一标识保护区的字段。

几何类型上要特别注意「多部件面」（MultiPolygon）。一个保护区可能由不相连的几块组成，比如某些湿地保护区包含多个湖泊。如果后续要做渔网分割或者按行政区统计，多部件面会导致面积重复计算或空间连接结果异常。判断方法很简单：在 QGIS 里看要素数量，再和属性表记录数对比，如果要素数少于记录数，说明存在多部件合并。

坐标系是第三个关键点。国内公开的保护区 shp 常见三种情况：WGS84 地理坐标系（EPSG:4326）、CGCS2000 地理坐标系（EPSG:4490）、以及各种高斯克吕格投影坐标系（如 CGCS2000 3 度带）。地理坐标系单位是度，直接算面积会得到平方度这种没有物理意义的数字；投影坐标系单位是米，才能算真实面积。很多人做面积统计时忘了这一步，结果出来的数字自己都不敢信。

2.2 用 Python 快速体检一份 shp 的完整脚本

下面这段代码用 geopandas 做一次完整体检，输出坐标系、几何类型、字段列表、面积范围和空几何数量。geopandas 的安装建议用 conda，避免 fiona 和 GDAL 的版本冲突。

import geopandas as gpd import pandas as pd # 读取 shp，注意路径不要有中文和空格 shp_path = "china_nature_reserves.shp" gdf = gpd.read_file(shp_path) # 1. 基本信息 print("要素数量:", len(gdf)) print("坐标系:", gdf.crs) print("几何类型:", gdf.geom_type.unique()) print("字段列表:", list(gdf.columns)) # 2. 检查空几何和无效几何 empty_count = gdf.geometry.is_empty.sum() invalid_count = (~gdf.geometry.is_valid).sum() print("空几何数量:", empty_count) print("无效几何数量:", invalid_count) # 3. 如果是地理坐标系，先粗算面积（仅用于对比，不精确） if gdf.crs and gdf.crs.is_geographic: # 用等面积投影估算，避免直接算平方度 gdf_proj = gdf.to_crs("EPSG:6933") # 全球等面积圆柱投影 gdf_proj["area_km2"] = gdf_proj.geometry.area / 1e6 print("面积范围(km2):", gdf_proj["area_km2"].min(), "-", gdf_proj["area_km2"].max()) print("面积总和(km2):", gdf_proj["area_km2"].sum()) # 4. 查看前几行属性 print(gdf.drop(columns="geometry").head())

这段代码的逻辑是：先读进来，把「有多少要素、什么坐标系、什么几何类型」这三件事确认清楚；然后检查空几何和无效几何，这两类问题在后续叠加分析时会导致崩溃或结果缺失；最后用 EPSG:6933 这个全球等面积投影做一次面积粗算，目的是快速判断数据量级是否合理，比如一个国家级保护区面积只有 0.1 平方公里，那大概率是坐标系或者单位出了问题。

参数说明：gdf.crs返回的是 pyproj 的 CRS 对象，如果显示 None，说明 shp 缺少 .prj 文件，需要手动指定坐标系。is_valid检查的是几何自相交、环方向错误等问题，无效几何可以用gdf.geometry.buffer(0)修复，但 buffer(0) 对复杂多部件面可能改变几何形状，修复后要重新检查面积变化。EPSG:6933 适合全球尺度，中国区域更推荐用 EPSG:4479（CGCS2000 大地坐标系）或按经度带选高斯克吕格投影。

2.3 坐标系不统一时怎么选投影

做全国尺度分析，我一般用 Albers 等面积投影，参数设中央经线 105°E，两条标准纬线 25°N 和 47°N，这是国内做全国土地利用、生态评估的常用投影。在 QGIS 里可以自定义 CRS，在 Python 里用 pyproj 拼字符串：

from pyproj import CRS albers_crs = CRS.from_proj4( "+proj=aea +lat_1=25 +lat_2=47 +lat_0=0 +lon_0=105 " "+x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs" ) gdf_albers = gdf.to_crs(albers_crs)

选 Albers 的理由是等面积特性，算出来的面积和真实面积偏差最小。如果做省级尺度，比如只分析四川或浙江的保护区，用对应经度带的高斯克吕格投影更合适，中央经线按 3 度带或 6 度带选。注意转换后要重新算面积，不要拿地理坐标系的面积字段直接用。

3. 从 shp 到可用分析图层：裁剪、合并与渔网分割的实操链路

3.1 按行政区裁剪保护区并统计各省面积

拿到全国保护区 shp 后，最常见的需求是按省统计保护区数量和面积。这里需要一份省级行政区划 shp 做空间连接或裁剪。我一般用空间连接（spatial join）而不是裁剪，因为裁剪会把跨省保护区切成两半，统计时容易重复计数。

# 读取省级行政区划 province = gpd.read_file("china_province.shp") # 统一坐标系到 Albers province_albers = province.to_crs(albers_crs) gdf_albers = gdf.to_crs(albers_crs) # 空间连接：每个保护区落在哪个省 joined = gpd.sjoin(gdf_albers, province_albers, how="left", predicate="intersects") # 按省统计数量和面积 stat = joined.groupby("PROVINCE").agg( reserve_count=("NAME", "nunique"), total_area_km2=("area_km2", "sum") ).reset_index() print(stat.sort_values("total_area_km2", ascending=False))

逻辑说明：sjoin的predicate="intersects"表示只要保护区和省份有交集就关联，跨省保护区会出现在多个省的结果里，所以用nunique去重计数。面积统计用sum会把跨省保护区面积在每个省都算一遍，如果需要精确到省的面积，应该先做交集（overlay intersection）再算面积。参数上，how="left"保证没有匹配到省份的保护区也保留，方便排查数据问题。

3.2 渔网分割 shp：把保护区切成规则网格做空间采样

渔网分割是热词里出现频率很高的操作，用途包括：把大范围保护区切成等大网格做抽样调查、生成训练样本、做空间异质性分析。QGIS 里有「创建网格」工具，Python 里用 shapely 和 geopandas 也能做。

import numpy as np from shapely.geometry import box def create_fishnet(gdf, cell_size_km=50): """在 gdf 范围内生成规则渔网""" minx, miny, maxx, maxy = gdf.total_bounds cell_size = cell_size_km * 1000 # 转米 cols = int(np.ceil((maxx - minx) / cell_size)) rows = int(np.ceil((maxy - miny) / cell_size)) cells = [] for i in range(cols): for j in range(rows): x0 = minx + i * cell_size y0 = miny + j * cell_size cells.append(box(x0, y0, x0 + cell_size, y0 + cell_size)) fishnet = gpd.GeoDataFrame({"cell_id": range(len(cells))}, geometry=cells, crs=gdf.crs) return fishnet # 生成 50km 渔网 fishnet = create_fishnet(gdf_albers, cell_size_km=50) # 用渔网裁剪保护区，得到每个网格内的保护区片段 clipped = gpd.overlay(gdf_albers, fishnet, how="intersection") clipped["clip_area_km2"] = clipped.geometry.area / 1e6 # 统计每个网格覆盖的保护区面积 grid_stat = clipped.groupby("cell_id")["clip_area_km2"].sum().reset_index() print(grid_stat.describe())

逻辑说明：total_bounds拿到整个图层的范围，按 50 公里步长生成正方形网格。overlay的how="intersection"做交集，把保护区切成网格内的碎片。clip_area_km2是每个碎片的面积，按cell_id汇总就得到每个网格覆盖的保护区总面积。

参数说明：cell_size_km根据分析尺度选，全国尺度一般 50 到 100 公里，省级尺度 10 到 20 公里，县级 1 到 5 公里。网格太小会导致碎片数量爆炸，内存扛不住；网格太大则失去空间采样意义。如果保护区分布极不均匀，可以考虑用六边形网格代替正方形，减少方向偏差。

3.3 shp 转 txt 和导出 wkt 的两种写法

热词里 shp 转 txt、shp 导出 wkt 的需求很实际，前者用于给不支持 shp 的系统喂数据，后者用于数据库存储或跨平台传输。转 txt 一般导出属性表加质心坐标，导出 wkt 则是把几何转成文本格式。

# 方式一：导出属性表 + 质心坐标到 txt gdf_albers["centroid_x"] = gdf_albers.geometry.centroid.x gdf_albers["centroid_y"] = gdf_albers.geometry.centroid.y gdf_albers.drop(columns="geometry").to_csv("reserves_centroid.txt", sep="\t", index=False, encoding="utf-8") # 方式二：导出 wkt 字符串 gdf_albers["wkt"] = gdf_albers.geometry.apply(lambda geom: geom.wkt) gdf_albers[["NAME", "wkt"]].to_csv("reserves_wkt.csv", index=False, encoding="utf-8")

逻辑说明：质心坐标适合做点状展示或空间插值，但注意质心不一定落在保护区内部，对于形状极不规则的保护区，质心可能跑到边界外。wkt 导出保留了完整几何信息，可以导入 PostGIS 或用于 Web 地图。参数上，sep="\t"用制表符分隔，避免保护区名称里的逗号干扰；encoding="utf-8"保证中文不乱码。如果 wkt 字符串太长，可以改用geom.simplify(0.01)做轻微简化，减小文件体积。

4. 避坑与排查：保护区 shp 处理中最容易翻车的五个地方

4.1 现象：面积统计结果比官方公布数字大好几倍

原因：坐标系是地理坐标系，直接对 geometry 调.area得到的是平方度，不是平方米。或者用了投影坐标系但单位是千米而非米，导致数量级错误。

解决：先确认gdf.crs的axis_info单位，地理坐标系必须转投影后再算面积。用gdf.crs.is_geographic判断，返回 True 就先to_crs到等面积投影。算完面积后和官方批复面积做对比，偏差超过 20% 就要检查投影参数和几何有效性。

4.2 现象：空间连接后同一个保护区在多个省重复出现

原因：跨省保护区在sjoin的intersects条件下会匹配多个省份，这是正常行为，但统计时没去重就会重复计数。

解决：计数用nunique而不是count，面积统计如果要求精确到省，先用overlay做交集再按省汇总。如果只是要一个「主要所属省」，可以按交集面积最大的省来归属。

4.3 现象：渔网分割后碎片数量异常多，程序卡死

原因：网格尺寸设得太小，或者保护区边界极其复杂（大量小多边形），交集运算产生海量碎片。

解决：先检查gdf.geom_type是否包含 MultiPolygon，用explode()拆成单部件面再分割。网格尺寸从大到小试，先 100 公里跑通流程再逐步缩小。如果碎片仍然过多，考虑用simplify对几何做轻微简化，容差取 0.001 度左右。

4.4 现象：shp 文件打开后属性表中文乱码

原因：shp 的 .dbf 文件默认编码是 GBK 或 Latin-1，而读取时用了 UTF-8。

解决：geopandas 读取时指定encoding="gbk"或encoding="latin1"试。如果还是乱码，用 QGIS 打开后重新导出为 UTF-8 编码的 GeoPackage 或 GeoJSON，再进 Python 处理。长期方案是统一转成 GeoPackage 格式，避免 shp 的编码和字段名长度限制。

4.5 现象：导出 wkt 后导入数据库报几何无效

原因：shp 里存在自相交、环方向错误或重复节点，wkt 原样导出后数据库的几何校验不通过。

解决：导出前用gdf.geometry.is_valid检查，对无效几何用buffer(0)或make_valid修复。修复后重新计算面积，确认几何没有发生大的形变。如果数据库是 PostGIS，可以用ST_MakeValid在入库时自动修复。

5. 进阶技巧：用保护区 shp 做空间自相关和热点分析

5.1 从渔网统计到空间权重矩阵

前面生成的渔网统计结果，可以直接拿来做空间自相关分析，判断保护区分布是聚集还是随机。核心是构建空间权重矩阵，然后算 Moran's I。Python 里用 libpysal 和 esda 两个库。

import libpysal from esda.moran import Moran # 基于渔网构建邻接权重 w = libpysal.weights.Queen.from_dataframe(fishnet) w.transform = "r" # 行标准化 # 把网格统计结果合并到渔网 fishnet = fishnet.merge(grid_stat, on="cell_id", how="left") fishnet["clip_area_km2"] = fishnet["clip_area_km2"].fillna(0) # 计算 Moran's I moran = Moran(fishnet["clip_area_km2"], w) print("Moran's I:", moran.I) print("p-value:", moran.p_sim)

逻辑说明：Queen.from_dataframe构建基于共享边或顶点的邻接关系，transform="r"做行标准化，让每个网格的邻居权重之和为 1。Moran计算全局自相关，I 大于 0 表示聚集，小于 0 表示离散，接近 0 表示随机。p_sim是置换检验的显著性，一般小于 0.05 才算显著。

参数说明：渔网尺寸会直接影响 Moran's I 的值，网格越大，空间聚集信号越容易被平滑掉。建议用多个尺寸（30km、50km、80km）分别算，看结果是否稳健。如果数据有大量零值网格（无保护区覆盖），可以考虑用局部 Moran's I（Local Moran）识别热点和冷点区域。

5.2 用局部 Moran 识别保护区热点县

全局 Moran's I 只给一个总体判断，实际规划中更关心「哪些地方是热点」。局部 Moran's I 可以输出每个网格的聚类类型：高-高聚集（热点）、低-低聚集（冷点）、高-低异常、低-高异常。

from esda.moran import Moran_Local lisa = Moran_Local(fishnet["clip_area_km2"], w) fishnet["lisa_cluster"] = lisa.q # 1=HH, 2=LH, 3=LL, 4=HL fishnet["lisa_p"] = lisa.p_sim # 筛选显著热点 hotspots = fishnet[(fishnet["lisa_cluster"] == 1) & (fishnet["lisa_p"] < 0.05)] print("显著热点网格数:", len(hotspots))

逻辑说明：lisa.q返回聚类类型编码，1 是高-高聚集，3 是低-低聚集。结合lisa_p < 0.05筛选统计显著的热点。这些热点网格对应的区域就是保护区分布最密集、最值得优先关注的地方。

参数说明：p_sim默认做 999 次置换，网格数量多时可以降到 499 次加快速度。如果热点网格连成大片，说明保护区分布有明显的空间聚集特征，可能和地形、气候带有关，可以叠加 DEM 或气候数据做进一步归因。

5.3 一个我常犯的错误

早期做保护区分析时，我习惯直接拿地理坐标系的 shp 算面积，结果一个国家级保护区算出 0.5 平方公里，差点当成数据错误报上去。后来养成习惯：拿到任何 shp，先跑一遍坐标系检查，再算面积，最后和官方数字对一遍。这个顺序不能省，省一步后面全是坑。另外渔网分割的网格尺寸不要一次定死，先用大网格跑通流程，确认代码没问题再缩小，不然调试时等交集运算能等到怀疑人生。希望帮到你。

本文还有配套的精品资源，点击获取