简介这份资源提供中国世界遗产地空间分布的Shapefile矢量数据面向地理信息系统学习者、空间分析研究者及需要中国遗产地底图数据的开发者可用于地图可视化、空间分布规律探究及与人口、经济等数据的叠加分析。压缩包共8个文件约21KB包含shp几何数据、dbf属性表、shx索引、prj坐标系统定义以及cpg编码说明、sbn与sbx空间索引和xml元数据等配套文件构成完整可用的GIS数据集。已有105人学习下载适合作为Python地理处理的练手素材。读者可借助geopandas、fiona等库直接读取快速绘制遗产地分布图统计各类遗产数量计算遗产地间空间距离并进一步结合其他数据源开展空间关联分析从而理解中国世界遗产的地理特征与分布规律。1. 拿到 shp 格式的中国世界遗产地空间分布数据先别急着打开 QGIS你手上如果有一个名为「中国世界遗产地空间分布」的压缩包解压后看到 .shp、.shx、.dbf、.prj 这一套文件说明你拿到的是一份标准的矢量点数据。它记录的是每一处世界遗产地的地理位置、名称、类别、列入年份等属性能直接用来做空间分布格局分析、核密度估计、缓冲区分析或者叠加到行政区划、地形、交通网络上做关联研究。适合谁做人文地理、遗产保护、旅游规划、城乡规划的研究生和从业者以及需要快速拿到一份可用点位数据做可视化或建模的人。但我要先泼一盆冷水shp 不是一个文件是一组文件少一个都打不开而且中文属性字段的编码问题几乎一定会让你翻车一次。这一章先把这份数据到底是什么、能干什么、坑在哪讲清楚后面再一步步带你跑通从解压到出图的全流程。2. 先搞懂 shp 到底是什么一组文件、三种编码、一个坐标系2.1 shp 不是单个文件是一套必须成套出现的文件组很多人第一次拿到 shp 数据看到压缩包里一堆同名不同后缀的文件就懵了。我一般会先列一下目录确认关键文件是否齐全。最小可用的组合是三个.shp 存几何形状.shx 存索引.dbf 存属性表。少了 .shx很多软件直接报错少了 .dbf你只能看到图形看不到任何名称和年份。如果要做投影转换或面积量算.prj 必须存在它定义了坐标系。还有一个 .cpg 文件专门声明 .dbf 的字符编码中文数据里这个文件极其关键没有它ArcGIS 和 QGIS 打开后属性表里的中文大概率变成乱码。# 列出解压后目录确认文件组是否齐全 ls -lh 中国世界遗产地空间分布/ # 典型输出应包含 # 中国世界遗产地空间分布.shp # 中国世界遗产地空间分布.shx # 中国世界遗产地空间分布.dbf # 中国世界遗产地空间分布.prj # 中国世界遗产地空间分布.cpg上面这条命令只是确认文件存在。如果发现缺 .prj你得去问数据提供方要坐标系定义或者根据数据范围自己判断——中国范围内的经纬度数据通常是 WGS84 地理坐标系EPSG:4326但也不绝对有些数据用的是 CGCS2000。缺 .cpg 的话可以自己新建一个同名 .cpg 文件里面写 UTF-8 或 GBK试一次就知道哪个对。2.2 坐标系和编码两个最容易让中文遗产地名变乱码的地方坐标系决定了你量算的距离和面积对不对。如果 .prj 里写的是 GCS_WGS_1984那你的单位是度直接算面积会得到平方度这种没有物理意义的数字。正确做法是先投影到适合中国的投影坐标系比如 Albers 等面积投影中央经线取 105°E双标准纬线取 25°N 和 47°N。这样量出来的面积单位才是平方米或平方公里。编码问题更隐蔽。.dbf 是 dBASE 格式它自己不声明编码靠 .cpg 告诉软件用 UTF-8 还是 GBK。我遇到过太多次QGIS 里打开属性表遗产名称全是问号或方块原因就是 .cpg 写的是 UTF-8但实际数据是 GBK 编码。解决办法有两个一是用 Python 的 geopandas 读取时显式指定 encodinggbk二是用 QGIS 的「重新编码」功能另存一份。下面这段代码演示如何安全读取并检查编码。import geopandas as gpd # 先尝试用 UTF-8 读取如果中文乱码再换 GBK try: gdf gpd.read_file(中国世界遗产地空间分布/中国世界遗产地空间分布.shp, encodingutf-8) # 检查名称字段是否包含乱码字符 sample gdf.iloc[0][名称] if 名称 in gdf.columns else gdf.iloc[0, 1] print(UTF-8 读取样例:, sample) except Exception as e: print(UTF-8 失败改用 GBK:, e) gdf gpd.read_file(中国世界遗产地空间分布/中国世界遗产地空间分布.shp, encodinggbk) print(GBK 读取样例:, gdf.iloc[0, 1]) # 查看坐标系和字段 print(坐标系:, gdf.crs) print(字段列表:, gdf.columns.tolist()) print(记录数:, len(gdf))这段代码的逻辑是先试 UTF-8失败或乱码再换 GBK。参数 encoding 直接传给底层读取库。注意 geopandas 读取 shp 时如果 .cpg 存在它会优先用 .cpg 里的编码所以有时候你指定了 encoding 也不生效这时候要么删掉 .cpg要么用 ogr2ogr 转换。字段列表里通常会有「名称」「类别」「列入年份」「所在省份」这几列具体以你拿到的数据为准不要假设字段名完全一致。2.3 属性表里有什么决定你后面能做什么分析一份典型的中国世界遗产地空间分布数据属性表至少包含名称、类别文化遗产/自然遗产/双重遗产、列入年份、所在省份。有些版本还会包含面积、缓冲区半径、是否跨境等字段。这些字段决定了你能做哪些分析有列入年份就能做时间序列的空间扩散分析有类别就能分类做核密度有省份就能按省汇总做统计图。如果字段缺失你得从其他来源补或者放弃对应分析。我一般会先用 pandas 把属性表导出来看一眼分布确认没有空值或异常值再往下做。import pandas as pd # 把属性表转成 DataFrame快速看分布 df pd.DataFrame(gdf.drop(columnsgeometry)) print(df.head()) print(df[类别].value_counts()) print(df[列入年份].describe())这里 drop(columnsgeometry) 是为了把几何列去掉只留属性。value_counts 看类别分布describe 看年份的统计量。如果发现年份有 0 或 9999 这种异常值说明数据有缺失标记后续分析要剔除。3. 用 Python 跑通从读取到出图的最小闭环3.1 环境准备geopandas 装不上怎么办geopandas 依赖 GDAL、Fiona、pyproj 这几个 C 库直接 pip install geopandas 在 Windows 上经常编译失败。我一般推荐用 conda 装省去编译麻烦。如果只能用 pip就先装预编译的 wheel 包。下面给出 conda 和 pip 两种方式按你的环境选一种。# 方式一conda推荐依赖自动解决 conda create -n heritage python3.10 conda activate heritage conda install -c conda-forge geopandas matplotlib contextily # 方式二pip需要先装 GDAL wheel pip install gdal3.6.2 --find-linkshttps://download.lfd.uci.edu/~gohlke/pythonlibs/ pip install geopandas matplotlib contextilyconda 方式最稳因为 conda-forge 频道里的 geopandas 已经把 GDAL 打包好了。pip 方式里那个 find-links 是预编译 wheel 的来源版本号要和你 Python 版本匹配3.10 就找 cp310 的包。装完后在 Python 里 import geopandas 不报错就说明环境通了。3.2 读取、投影转换、核密度分析一条龙拿到数据后我通常先做三件事读进来、转投影、做核密度。核密度能直观看出遗产地在空间上是不是聚集聚集在哪里。下面这段代码把这三步串起来最后出一张核密度图。import geopandas as gpd import matplotlib.pyplot as plt from shapely.geometry import Point import numpy as np # 1. 读取数据处理编码 gdf gpd.read_file(中国世界遗产地空间分布/中国世界遗产地空间分布.shp, encodinggbk) # 2. 如果坐标系是地理坐标系先投影到 Albers 等面积投影 if gdf.crs and gdf.crs.is_geographic: gdf gdf.to_crs(projaea lat_125 lat_247 lat_00 lon_0105 datumWGS84 unitsm) # 3. 提取点坐标做核密度 coords np.array([(geom.x, geom.y) for geom in gdf.geometry]) from scipy.stats import gaussian_kde kde gaussian_kde(coords.T) # 生成网格 xmin, ymin, xmax, ymax gdf.total_bounds xx, yy np.mgrid[xmin:xmax:200j, ymin:ymax:200j] positions np.vstack([xx.ravel(), yy.ravel()]) density kde(positions).reshape(xx.shape) # 4. 出图 fig, ax plt.subplots(figsize(10, 8)) ax.imshow(np.rot90(density), cmapYlOrRd, extent[xmin, xmax, ymin, ymax]) gdf.plot(axax, markersize5, colorblack, alpha0.6) ax.set_title(中国世界遗产地核密度分布) plt.savefig(heritage_kde.png, dpi300, bbox_inchestight) plt.show()这段代码的关键点to_crs 里的投影参数是 Albers 等面积投影适合中国全境lat_1 和 lat_2 是双标准纬线lon_0 是中央经线。gaussian_kde 的带宽默认用 Scott 规则如果觉得太平滑或太尖锐可以手动传 bw_method0.2 调整。total_bounds 拿到的是投影后的范围单位是米。imshow 里 np.rot90 是因为 mgrid 生成的数组方向和图的方向差 90 度。最后保存成 300 dpi 的 PNG直接能放进论文或报告。3.3 叠加行政区划和底图让分布图能讲故事光有点和核密度还不够评审或读者需要看到遗产地和省界、地形、水系的关系。叠加行政区划最简单的方式是读一份省界 shp然后 plot 上去。如果没有省界数据可以用 contextily 加在线底图但要注意在线底图需要网络且坐标系必须是 Web MercatorEPSG:3857。下面演示叠加省界的做法。# 假设你有一份省界 shp坐标系也是 WGS84 province gpd.read_file(省界数据/provinces.shp, encodingutf-8) province province.to_crs(gdf.crs) # 统一到同一投影 fig, ax plt.subplots(figsize(12, 10)) province.plot(axax, facecolornone, edgecolorgray, linewidth0.5) gdf.plot(axax, markersize8, colorred, alpha0.7) ax.set_title(中国世界遗产地与省界叠加) plt.savefig(heritage_province.png, dpi300, bbox_inchestight)province.to_crs(gdf.crs) 这一步必须做否则两个图层坐标系不一致叠上去位置会偏到太平洋。facecolornone 让省界只显示边框不填充颜色避免盖住下面的点。markersize 和 alpha 根据点的密度调点太密就调小 markersize 或加大 alpha 透明度。4. 避坑指南中文乱码、坐标系错位、字段丢失的排查清单4.1 中文属性全是乱码改了 .cpg 也没用现象QGIS 或 ArcGIS 打开属性表遗产名称显示为「????」或「锟斤拷」。原因.dbf 实际编码和 .cpg 声明不一致或者软件读取时忽略了 .cpg。解决用 Python 显式指定 encoding 读取确认正确编码后用 ogr2ogr 重新导出并强制写入 .cpg。命令如下ogr2ogr -f ESRI Shapefile output.shp input.shp -lco ENCODINGUTF-8-lco ENCODINGUTF-8 会同时写 .cpg 文件确保后续打开不乱码。如果原始数据是 GBK先读进来再导出为 UTF-8不要直接改 .cpg 内容因为 .dbf 里的字节没变改声明没用。4.2 点位置整体偏移跑到海里去了现象叠加省界后发现遗产点整体偏移几百米到几公里。原因两个图层坐标系不一致或者其中一个坐标系定义错误。解决先用 gdf.crs 打印坐标系确认是否一致。如果不一致用 to_crs 统一。如果 .prj 文件缺失或错误需要手动指定 crs 后再转换。常见错误是把 WGS84 的数据当成 CGCS2000 用两者在中国范围内差异不大但严格来说不能混用。4.3 字段名变成乱码或截断现象属性表里字段名显示为「Name_1」「Field_2」这种默认名。原因.dbf 字段名长度限制是 10 个字符中文字段名在导出时可能被截断或转码失败。解决在 Python 里读取后重命名字段或者用 ogr2ogr 导出时用 -sql 指定字段别名。重命名示例gdf gdf.rename(columns{名称: name, 类别: category, 列入年份: year})重命名后再做分析避免后续代码里引用中文字段名出问题。4.4 核密度图一片空白或全黑现象核密度图要么没有颜色要么全黑。原因网格分辨率太低或太高或者密度值范围极端。解决调整 mgrid 的步长200j 改成 500j 或 100j或者对密度值做对数变换。另外检查 coords 是否为空如果 gdf 没有几何数据coords 就是空数组kde 会报错。4.5 保存的图片边缘被裁掉现象savefig 出来的图标题或图例被切掉。原因bbox_inchestight 有时候会裁过头。解决改用 bbox_inchesstandard或者手动调 subplots_adjust。我一般用 plt.tight_layout() 配合 bbox_inchestight大部分情况够用。5. 进阶技巧用空间自相关验证遗产地是否真的聚集5.1 全局莫兰指数一个数字判断聚集还是分散核密度图是视觉判断要定量验证聚集性用全局莫兰指数Global Morans I。它输出一个介于 -1 到 1 之间的值大于 0 表示聚集小于 0 表示分散接近 0 表示随机。下面用 libpysal 和 esda 计算。import libpysal from esda.moran import Moran import numpy as np # 基于点坐标构建空间权重矩阵用 K 近邻 coords np.array([(geom.x, geom.y) for geom in gdf.geometry]) w libpysal.weights.KNN.from_array(coords, k8) w.transform r # 行标准化 # 用某个数值字段做自相关比如列入年份 y gdf[列入年份].values moran Moran(y, w) print(fMorans I: {moran.I:.4f}, p-value: {moran.p_sim:.4f})KNN 的 k8 是经验值点少可以调小点多可以调大。w.transformr 做行标准化让每个点的邻居权重和为 1。moran.p_sim 是置换检验的伪 p 值小于 0.05 说明聚集显著。如果 I 为正且显著说明列入年份相近的遗产地在空间上趋于聚集可能反映某种区域申报规律。5.2 局部莫兰指数找出热点和冷点区域全局指数只给一个总数局部莫兰指数Local Morans I能告诉你哪些具体位置是热点高-高聚集、冷点低-低聚集或异常值。用 esda 的 Moran_Local 计算然后画 LISA 聚类图。from esda.moran import Moran_Local import matplotlib.pyplot as plt lisa Moran_Local(y, w) fig, ax plt.subplots(figsize(10, 8)) gdf.plot(axax, columnlisa.q, categoricalTrue, legendTrue, markersize20) ax.set_title(局部莫兰指数聚类图q1 高-高q3 低-低) plt.savefig(lisa_cluster.png, dpi300, bbox_inchestight)lisa.q 是聚类类型1 代表高-高2 代表低-高3 代表低-低4 代表高-低。画图时用 categoricalTrue 让不同类别显示不同颜色。这张图能直接指出哪些省份或区域是遗产地申报的热点区哪些是冷点区比单纯核密度图更有解释力。5.3 我踩过的坑和固定习惯做这份数据的过程中我最大的教训是不要相信任何 shp 文件自带的 .cpg 和 .prj一定要自己用代码读一遍、打印坐标系、检查字段编码。我现在的固定习惯是拿到任何 shp 先跑三行代码print(gdf.crs)、print(gdf.columns)、print(gdf.head(3))。这三行能提前暴露 90% 的问题。另外核密度分析的带宽不要用默认值手动调两三次对比效果再定。空间自相关的权重矩阵KNN 的 k 值也要试k4、6、8 各跑一次看结果稳不稳定。这些参数没有绝对标准但多试几次能让你对数据的空间结构更有感觉。希望帮到你。本文还有配套的精品资源点击获取