企业官网建设流程全解析

简介：这份中国自然保护区空间分布数据以shp矢量格式整理，面向从事GIS制图、生态地理研究与空间分析的用户，可用于全国及区域尺度的保护区分布制图、叠加分析与规划参考。资源按国家级与地方级两类保护区分别组织，每类均配套完整的shapefile组件，包括存储几何要素的shp、记录属性字段的dbf、定义坐标系统的prj，以及shx、sbn、sbx、cpg、xml等索引与元数据文件，便于在ArcGIS、QGIS等平台直接加载、符号化与空间查询。压缩包共16个文件，整体约7.95MB，体量轻便，适合课堂练习、论文插图与项目底图使用。目前已有948人学习下载，说明该数据在生态保护与地理信息相关场景中具有一定参考价值，能帮助读者省去自行数字化与配准的工作量，快速获得可直接参与分析的分级保护区图层。

1. 从一份全国自然保护区 shp 说起：谁在找它，找到后能干什么

如果你正在做生态红线评估、国土空间规划底图，或者要给某个区域做生物多样性廊道分析，大概率绕不开“自然保护区边界”这层数据。中国自然保护区空间分布数据 shp 就是干这个的——它把全国各级自然保护区的位置和范围，以矢量面要素的形式存成 ESRI Shapefile，能在 ArcGIS、QGIS、PostGIS 里直接打开、叠加、裁剪。我第一次拿到这类数据是为了给一个风电场选址做生态敏感区避让，当时手上只有 DEM 和行政区划，缺保护区边界，整个分析卡住。这份 shp 解决的就是“保护区在哪、边界长什么样、和我的项目区有没有重叠”这三个问题。适合做国土、生态、遥感、环境咨询的从业者，也适合刚接触 GIS 矢量数据、想拿真实数据练手的新手。它不是瓦片底图，不是遥感影像，是实打实的矢量面，能参与空间运算，这才是它值钱的地方。

2. 拿到 shp 先别急着打开：坐标系、字段和拓扑三件事

2.1 坐标系不确认，面积和距离全是错的

Shapefile 本身可以带 .prj 文件，也可以不带。很多流传的保护区 shp 是地理坐标系（GCS_WGS_1984 或 GCS_China_Geodetic_Coordinate_System_2000），单位是度，不是米。你直接用它算面积，得到的是“平方度”，没有任何物理意义。常见做法是先确认 .prj 内容，再决定要不要投影。我一般会在 QGIS 里右键图层属性看 CRS，或者用 ogrinfo 命令行查。

# 查看 shp 的坐标系和字段结构，ogrinfo 来自 GDAL ogrinfo -so -al china_nature_reserves.shp # 输出里重点看两行： # Geometry: Polygon # Feature Count: xxx # 以及 Coordinate System 后面的 EPSG 编号

逻辑说明：-so只输出摘要不输出每个要素，-al列出所有图层。参数上，如果 Coordinate System 显示GCS_WGS_1984且没有 EPSG 投影代号，说明它是地理坐标。要算面积，得先转成投影坐标，比如 Albers 等积投影或 UTM 对应带号。我一般会统一转成EPSG:4526（CGCS2000 / 3-degree Gauss-Kruger zone 39）或EPSG:3857做可视化，但算面积一定用等积投影，别用 Web Mercator，它的面积畸变在高纬度能到几倍。

2.2 字段里藏着保护区级别和名称，别只拿几何

保护区 shp 的属性表通常有NAME、LEVEL、TYPE、PROVINCE这类字段，但不同来源字段名不一样。有的用中文，有的用拼音缩写，有的干脆只有Id和name。打开后先看字段，再决定怎么筛选。比如你要只保留国家级，就得找LEVEL或GRADE字段，值可能是“国家级”“National”或数字编码。

# 用 geopandas 读取并查看字段，比在 ArcGIS 里点属性表快 import geopandas as gpd gdf = gpd.read_file("china_nature_reserves.shp") print(gdf.columns.tolist()) print(gdf[['NAME', 'LEVEL']].head(10) if 'LEVEL' in gdf.columns else gdf.head(10)) print("要素总数:", len(gdf)) print("几何类型:", gdf.geom_type.unique())

逻辑说明：gpd.read_file直接读 shp，返回 GeoDataFrame。columns.tolist()看字段名，head(10)抽样看值。如果LEVEL字段不存在，就用gdf.columns里找含“级”“level”“grade”的列。参数上，len(gdf)是要素数，geom_type.unique()确认是不是全是 Polygon，如果混了 MultiPolygon 或 Point，后面做叠加会出问题。

2.3 拓扑检查：自相交和重叠面会让后续分析翻车

保护区边界数据常见的问题是面自相交、相邻保护区重叠、或者有极小的碎片多边形。这些在可视化时看不出来，一做相交分析或面积统计就报错。我一般会先跑一遍有效性检查。

# 检查几何有效性，修复自相交 from shapely.validation import explain_validity invalid = gdf[~gdf.is_valid] print("无效几何数量:", len(invalid)) if len(invalid) > 0: print(invalid.geometry.apply(explain_validity).head()) # 修复：用 buffer(0) 是常见做法，但会轻微改变边界 gdf_fixed = gdf.copy() gdf_fixed['geometry'] = gdf_fixed.geometry.buffer(0) print("修复后无效数量:", len(gdf_fixed[~gdf_fixed.is_valid]))

逻辑说明：is_valid返回布尔序列，explain_validity告诉你为什么无效，比如“Self-intersection”。buffer(0)是经典修复手段，但它会把自相交的环重新组织，可能产生微小位移。如果精度要求高，用make_valid或 ArcGIS 的 Repair Geometry 工具。参数上，buffer(0)的距离是 0，只做拓扑修复不做缓冲扩张。修复后一定要再统计一次面积，和修复前对比，差异超过 1% 就要警惕。

3. 把保护区 shp 用起来：叠加、裁剪、导出和格式转换

3.1 和项目区做相交：判断是否涉及保护区

这是最常见的需求。你有一个项目边界 shp，想知道它和保护区有没有重叠、重叠多少。核心操作是空间相交。

import geopandas as gpd # 读取项目区和保护区 project = gpd.read_file("project_area.shp") reserves = gpd.read_file("china_nature_reserves.shp") # 统一坐标系，以项目区为准 if project.crs != reserves.crs: reserves = reserves.to_crs(project.crs) # 相交分析 intersect = gpd.overlay(project, reserves, how='intersection') print("重叠要素数:", len(intersect)) if len(intersect) > 0: intersect['overlap_area'] = intersect.geometry.area print(intersect[['NAME', 'overlap_area']].sort_values('overlap_area', ascending=False))

逻辑说明：to_crs统一坐标系，避免“看起来重叠但实际错位”的玄学问题。gpd.overlay的how='intersection'保留两者重叠部分。overlap_area是投影坐标系下的面积，单位取决于 CRS，如果是米制就是平方米。参数上，如果项目区是点数据，用sjoin更合适；如果是面，用overlay。注意overlay对无效几何很敏感，前面修复过的数据在这里才稳。

3.2 按级别或省份筛选后导出

不是每次都要全国数据。做省级规划时，只保留目标省份的保护区，文件小、加载快。

# 按省份字段筛选，字段名可能是 PROVINCE 或省 province_field = 'PROVINCE' if 'PROVINCE' in reserves.columns else '省' target = reserves[reserves[province_field].str.contains('四川', na=False)] # 导出为新的 shp，注意编码 target.to_file("sichuan_reserves.shp", encoding='utf-8') print("导出要素数:", len(target))

逻辑说明：str.contains('四川')做模糊匹配，na=False避免空值报错。to_file的encoding='utf-8'很关键，Shapefile 的 .dbf 默认编码在中文环境容易乱码，加上这个参数能减少翻车概率。参数上，如果字段值是“四川省”而不是“四川”，contains依然能命中。导出后建议用 ogrinfo 再确认一次要素数和坐标系。

3.3 shp 转 GeoJSON、WKT 和 3D Tiles 的边界

热词里出现 dwg转shp、shp转3dtiles、shp转txt、shp导出wkt，说明大家有格式转换需求。这里说清楚：shp 转 GeoJSON 和 WKT 是常规操作，转 3D Tiles 是另一条技术路线，需要先有高度信息或建筑体块，单纯保护区面数据转 3D Tiles 只能做贴地展示，没有立体意义。

# shp 转 GeoJSON reserves.to_file("reserves.geojson", driver='GeoJSON') # 导出每个要素的 WKT，用于数据库或文本交换 reserves['wkt'] = reserves.geometry.apply(lambda g: g.wkt) reserves[['NAME', 'wkt']].to_csv("reserves_wkt.csv", index=False, encoding='utf-8-sig') # shp 转 txt：通常是把坐标点导出，面数据转 txt 意义有限 # 如果是点数据，可以这样 # points = gpd.read_file("points.shp") # points['x'] = points.geometry.x # points['y'] = points.geometry.y # points[['x', 'y']].to_csv("points.txt", sep='\t', index=False)

逻辑说明：driver='GeoJSON'指定输出格式，GeoJSON 默认 WGS84，如果原数据是投影坐标，to_file不会自动转，需要先to_crs(4326)。WKT 导出用g.wkt，适合入库或传给其他系统。encoding='utf-8-sig'让 Excel 打开 CSV 不乱码。参数上，sep='\t'是制表符分隔，有些老系统只认这种 txt。注意：shp 转 3D Tiles 不在这个脚本范围，那需要 Cesium ion 或专门工具链，且需要高度字段。

4. 避坑与排查：保护区 shp 最常见的五个翻车现场

4.1 打开后一片空白，什么也看不见

现象：在 ArcGIS 或 QGIS 里加载 shp，图层列表有，但地图窗口空白。原因：坐标系未定义或定义错误，软件把坐标当成经纬度但实际是投影坐标，或者反过来，导致视图缩放到错误位置。解决：先ogrinfo看 CRS，如果缺失，用gdf.set_crs(epsg=4326, allow_override=True)强制指定，再to_crs到正确坐标系。如果 CRS 明显不对，比如数据在 China 但显示在非洲，就是定义错了，重新指定。

4.2 面积算出来是 0 或者负数

现象：用geometry.area算面积，结果全是 0 或极小值。原因：数据是地理坐标系，单位是度，面积数值极小；或者几何无效导致面积为 0。解决：先to_crs到等积投影，比如 Albers 或 UTM，再算面积。同时检查is_valid，无效几何先修复。我一般会固定用EPSG:4526或EPSG:3857做可视化，但面积统计一定用等积投影，这是血泪经验。

4.3 中文属性乱码

现象：属性表里保护区名称显示成“?????”或乱码。原因：Shapefile 的 .dbf 文件编码和软件读取编码不一致，常见于 GBK 和 UTF-8 混用。解决：用 geopandas 读取时指定encoding='gbk'或encoding='utf-8'试，导出时统一用encoding='utf-8'。如果已经乱码，用gdf['NAME'] = gdf['NAME'].str.encode('latin1').str.decode('gbk')尝试修复，但不保证成功。预防办法是拿到数据先看 .cpg 文件，里面写了编码。

4.4 叠加分析报拓扑异常

现象：gpd.overlay或 ArcGIS Intersect 报“TopologyException”或“invalid geometry”。原因：面自相交、重复节点、零面积碎片。解决：先buffer(0)修复，或者用make_valid。如果修复后还报错，把数据拆成小块逐块排查，定位到具体要素后手动修。常见做法是导出无效要素，在 QGIS 里用“顶点工具”手动改，虽然慢但准。

4.5 筛选省份时匹配不到

现象：用str.contains('四川')筛选，结果为空。原因：字段值可能是“四川省”“Sichuan”“SC”或带空格。解决：先print(gdf['PROVINCE'].unique())看实际值，再决定匹配规则。如果值是拼音，就用拼音匹配；如果有空格，先str.strip()。参数上，str.contains默认区分大小写，加case=False更稳。

5. 进阶：用渔网分割和空间连接做保护区覆盖率统计

5.1 渔网分割 shp：把保护区切成规则网格

热词里有“渔网分割shp”，这在做空间统计时很实用。比如你要统计每个 10km×10km 网格内的保护区面积占比，就得先建渔网，再和保护区相交。

import geopandas as gpd from shapely.geometry import box import numpy as np # 读取保护区并投影到米制 reserves = gpd.read_file("china_nature_reserves.shp").to_crs(epsg=4526) # 获取范围 minx, miny, maxx, maxy = reserves.total_bounds # 创建 10km 渔网 cell_size = 10000 # 米 cols = np.arange(minx, maxx, cell_size) rows = np.arange(miny, maxy, cell_size) cells = [box(x, y, x + cell_size, y + cell_size) for x in cols for y in rows] grid = gpd.GeoDataFrame({'geometry': cells}, crs=reserves.crs) # 相交计算每个网格的保护区面积 intersect = gpd.overlay(grid, reserves, how='intersection') intersect['area'] = intersect.geometry.area coverage = intersect.groupby(intersect.index)['area'].sum().reset_index() grid['reserve_area'] = grid.index.map(coverage.set_index('index')['area']).fillna(0) grid['coverage_ratio'] = grid['reserve_area'] / (cell_size ** 2) # 导出结果 grid.to_file("reserve_coverage_grid.shp", encoding='utf-8') print("网格总数:", len(grid)) print("有保护区覆盖的网格数:", len(grid[grid['coverage_ratio'] > 0]))

逻辑说明：total_bounds拿到保护区范围，box生成每个网格。overlay相交后按网格索引分组求和，得到每个网格内的保护区面积。coverage_ratio是覆盖率，0 到 1 之间。参数上，cell_size按需改，做省级用 10km，做全国用 50km 更合适。注意grid.index.map那里，coverage的index列是原网格索引，映射回去后fillna(0)处理没有重叠的网格。这个脚本跑全国数据可能慢，建议先按省筛选。

5.2 空间连接：给每个网格打上保护区名称标签

如果你想知道每个网格涉及哪些保护区，用sjoin做空间连接。

# 空间连接，保留网格和保护区名称 joined = gpd.sjoin(grid, reserves[['NAME', 'geometry']], how='left', predicate='intersects') # 每个网格可能对应多个保护区，聚合名称 name_agg = joined.groupby(joined.index)['NAME'].apply(lambda x: '|'.join(x.dropna().unique())).reset_index() grid['reserve_names'] = grid.index.map(name_agg.set_index('index')['NAME']).fillna('无') print(grid[grid['reserve_names'] != '无'][['reserve_names', 'coverage_ratio']].head())

逻辑说明：sjoin的predicate='intersects'表示只要几何有接触就连接。how='left'保留所有网格，没有保护区的网格NAME为 NaN。聚合时用'|'.join把多个保护区名拼成一个字段，方便查看。参数上，predicate还可以用within或contains，取决于你想让网格“包含”保护区还是“相交”保护区。这一步做完，你就能回答“哪些网格涉及国家级保护区”这种问题了。

5.3 验证数据质量的三个习惯

从那以后我每次拿到新的 shp，都强制走一遍这三步：第一，ogrinfo看 CRS 和要素数；第二，is_valid查几何有效性；第三，抽样对比面积和官方公布值。比如某个保护区的面积如果和公开资料差 10% 以上，就要怀疑边界精度或坐标系问题。这套习惯帮我拦下过好几次“数据看起来能用但一算就错”的情况。希望帮到你。

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