尧图网页设计 
简介面向户外运动研究、城市规划与智慧城市等领域这份GIS数据集提供了广州市2020年徒步活动的完整轨迹记录。数据以标准Shapefile格式存储包含轨迹ID、徒步距离、运动速度、采集日期等核心属性可广泛应用于市民徒步行为分析、热门路线识别、运动设施需求评估等场景适合地理信息数据分析师、城市规划师及运动科学方向的研究者使用。压缩包共7个文件主文件为shp几何数据配套的dbf属性表保存了全部轨迹字段prj文件定义了坐标系统xml元数据提供了数据集说明sbn、sbx、shx索引文件则用于加速空间查询与读取包体约649.82MB。目前已有265人学习或下载。每个轨迹ID可区分不同徒步活动距离与速度属性为运动强度分析提供依据直接使用可省去自行采集GPS轨迹的繁琐流程快速开展徒步空间分布、区域运动资源配置等研究结合天气、人口密度等多源数据还能进行多因素深度挖掘为城市公共空间优化与公共健康决策提供数据支撑。1. 把一份轨迹shp当资产而不是图层收到“广州市户外运动轨迹-2020年徒步轨迹shp”这类数据时第一反应通常是拖进ArcMap看一眼但几秒钟后就会发现满屏都是密密麻麻的折线属性表里全是拼音缩写根本分不清哪条是白云山、哪条是火炉山。这份数据的价值也不在“能显示出来”而在回答“哪段路被走得最多”“哪些线路一年里重复率最高”“热度分布围绕什么山体展开”这类实际问题。shp本身不是轨迹专用格式它把带有时间语义的轨迹压成了二维线要素属性字段还经常混着起点终点所以直接分析会连续踩坐标系、拓扑和字段语义三个坑。下面这套处理流程专为轨迹型shp设计用GeoPandas做读取、清洗和统计用GDAL和QGIS做导出与批量出图不依赖ArcGIS也不需要企业级GIS许可适合户外赛事组织者、城市绿道评估和路网优化方向的工程师。2. 先建立坐标系与字段映射别急着画线轨迹shp最常见的问题不是几何错误而是坐标系声明缺失或字段单位不一致。打开要素前先问自己三个问题坐标是经纬度还是投影米每条线是一条完整的徒步记录还是多个线段叠加属性表里的时间字段能不能直接解析这三件事不解决后面的抽稀、测长、渔网统计全部会失真。2.1 用GeoPandas读取并核对轨迹shp的几何类型先用一段最小化代码把shp的关键信息打印出来这一步能避免后面80%的错误。import geopandas as gpd trails gpd.read_file(广州市户外运动轨迹-2020年徒步轨迹.shp) print(trails.crs) print(trails.geom_type.unique()) print(trails.columns.tolist()) print(trails.head(3))geom_type.unique()会给出LineString或MultiLineString两种典型结果。如果是后者后续length属性和坐标遍历都需要按多部件逐段处理。crs输出None意味着shp缺少.prj文件GeoPandas会按无投影坐标读入此时算出的长度是十进制度数而不是米。columns.tolist()用来确认字段名里是否有start_time、end_time、name、length如果没有时间字段后面的时序拆分就没法做。拿到字段名后还要看一眼坐标范围判断数据是否真的落在广州附近。print(trails.total_bounds)广州的经纬度范围大致是东经113.0至114.5北纬22.5至24.0。如果total_bounds的小数点前是六位数比如[380000, 2600000, 390000, 2620000]说明这份shp已经做过投影不能直接按WGS84处理。2.2 统一坐标系先声明EPSG:4326再投到米制轨迹采集设备输出的原生坐标通常是WGS84经纬度ESRI的shp组件会把投影信息写进.prj但单独压缩传递时经常漏掉这个文件。遇到crs为None的情况不要盲目使用to_crs而是先声明。if trails.crs is None: trails trails.set_crs(EPSG:4326) else: print(已有坐标系, trails.crs)set_crs和to_crs的区别一定要分清前者只是给坐标贴上坐标系标签不改变数值后者才是真正做投影变换。如果坐标范围显示是投影坐标我却把它声明成4326后面的所有长度都会错得离谱。稳妥的判断方式是看total_bounds的绝对量级经纬度在0到180之间米制在10万到上千万之间。测距和做渔网统计前需要转换到适合广州的投影坐标系。trails_m trails.to_crs(EPSG:32649) trails_m[length_m] trails_m.geometry.length print(trails_m[length_m].describe())EPSG:32649是WGS 84 / UTM zone 49N广州恰好落在中央经线117°W附近长度变形很小适合做距离计算和栅格统计。不同坐标系各有用处不是所有米制坐标都适合算长度。EPSG名称单位主要用途4326WGS 84度原始轨迹存储、KML/TXT交换、GPS设备输出32649WGS 84 / UTM zone 49N米长度测量、缓冲区、渔网统计、核密度分析3857Web墨卡托米在线底图叠加显示不适合面积和距离计算坐标转换后再看length_m的统计值如果发现中位数只有几米说明原始shp里的线已经被打断需要先按轨迹ID合并而不是做分析。2.3 时间字段解析把shp变成可按时序分析的数据轨迹shp属性表里常见时间字段有start_time、end_time、date。这些字段在Excel里看着正常读进GeoPandas后可能还是字符串分词需要显式转换。import pandas as pd trails[start_dt] pd.to_datetime( trails[start_time], errorscoerce ) trails[year] trails[start_dt].dt.year trails[month] trails[start_dt].dt.month print(trails[month].value_counts().sort_index())errorscoerce表示解析失败时置为NaT不会因为个别脏数据中断整个脚本。解析完成后马上统计month分布可以判断这条shp是否真的覆盖2020年全年还是只集中在下半年。如果字段原本是时间戳数字比如1609430400需要先除以86400再转日期不能直接丢给pd.to_datetime。时间字段解析完成才谈得上清理几何。很多轨迹shp在采集时没有做航点抽稀一条白云山环线可能包含上万个顶点几何修复的优先级反而高于抽稀。3. 清洗轨迹shp几何修复、抽稀和字段拆分轨迹shp的几何脏点通常有三类自相交导致is_validFalse大量冗余顶点拖慢绘制和计算同一条轨迹被GPS中断拆成多个要素。这三类问题不解决后续聚合统计的结果都会偏高或偏低。清洗的目标不是把几何改得多整齐而是让它满足后续“按线聚合、按米求和”的基本前提。3.1 修复无效几何与自相交线段用Shapely自带的make_valid可以一次性处理大部分几何问题。from shapely.validation import make_valid from shapely.ops import linemerge trails[geometry] trails.geometry.apply(make_valid) def to_single_line(geom): if geom.geom_type MultiLineString: merged linemerge(geom) return merged if merged.geom_type LineString else geom return geom trails[geometry] trails.geometry.apply(to_single_line)make_valid会把自相交的线拆成多个部件所以调用后要用linemerge把首尾相接的部件重新合并成一条完整线。linemerge只处理端点完全重合的情况如果GPS断点之间有几米空隙接口会保持MultiLineString这一步也算变相检测出轨迹中断。修复完成后检查无效几何比例print(无效比例, round((~trails.is_valid).mean(), 4))如果无效比例超过1%说明数据源质量较差后面统计网格时要把无效要素先过滤掉否则overlay可能报错。3.2 抽稀用simplify控制坐标点密度轨迹shp的顶点往往每隔一两米就有一个出图和计算都容易被数据量拖垮。道格拉斯-普克算法的目的不是简化轮廓而是去掉对几何形状影响小于阈值的点。trails_m[geom_simp] trails_m.geometry.simplify( tolerance2.0, preserve_topologyTrue ) trails_clean trails_m.set_geometry(geom_simp).drop(columns[geometry])tolerance单位与坐标系一致。前面我们用了EPSG:32649所以2.0代表2米。容差越大顶点越少但也越容易把急转弯拉直。tolerance适用场景0.5米保留栈道拐角、观景台停留点1.0米普通山地徒步轨迹2.0米广州山体绿道、城市步道5.0米公里级热度底图不适合作精确路线抽稀后记得重新计算长度因为简化会稍微缩短线长。若原始轨迹点非常密集抽稀前后总长度差异应小于1%差异过大说明容差设得太大。3.3 缺失字段拆分与分组输出清洗完几何后很多场景需要按月份或路线名称输出子集。常见的做法是按month字段分组逐组写出独立的shp文件。for month, group in trails_clean.groupby(month): group.to_file( ftrails_2020_{month:02d}.shp, encodingutf-8 ) print(month, len(group), group[length_m].sum())这里的encodingutf-8只影响.dbf属性表。如果后续用ArcGIS打开老版本读UTF-8中文可能乱码一般建议在ArcGIS里也设置代码页。实际上更稳的交付方案是保留字段名用拼音或英文把中文名称放到name字段值里。字段拆分还可以按轨迹距离过滤。比如只保留长度大于1公里的轨迹用来过滤掉GPS漂移产生的碎线。trails_clean trails_clean[trails_clean[length_m] 1000]轨迹清洗阶段不需要过多纠结“线是否光滑”重点是让每个要素都能代表一次独立徒步。下一步要做的是把这些线变成可以排序的密度指标。4. 用渔网和核密度估算徒步热点路线是否热门不能只靠“看起来线多”来判断。把轨迹shp切割到规则网格里统计每个网格内轨迹线的累计长度比目测叠加更客观。这一节先用渔网做分段统计再用核密度做连续热力面。4.1 在投影坐标系下建立200米渔网渔网大小直接影响统计粒度。广州山体徒步轨迹的路线宽度通常在1到3米网格设到100米会切碎同一条山脊线设到500米又分不清白云山和火炉山的边界。200米是折中值。from shapely.geometry import box xmin, ymin, xmax, ymax trails_clean.total_bounds cellsize 200 cols int((xmax - xmin) // cellsize) 1 rows int((ymax - ymin) // cellsize) 1 grid_polys [] for i in range(cols): for j in range(rows): x0 xmin i * cellsize y0 ymin j * cellsize grid_polys.append(box(x0, y0, x0 cellsize, y0 cellsize)) grid gpd.GeoDataFrame(geometrygrid_polys, crstrails_clean.crs) grid[grid_id] range(len(grid))这段代码根据轨迹总范围生成一个全覆盖规则网格每个格子都是独立的矩形要素。grid_id是后续聚合的关键索引不要省略。接下来做线和网格的相交计算。inter gpd.overlay(grid, trails_clean[[geometry]], howintersection) inter[seg_len] inter.geometry.length density ( inter.groupby(grid_id)[seg_len] .sum() .reset_index() ) grid_density grid.merge(density, ongrid_id, howleft) grid_density[seg_len] grid_density[seg_len].fillna(0)overlay会把落在多个网格里的轨迹线逐一分割seg_len是每段在网格内的实际长度。按grid_id汇总后grid_density就带上了每个网格的累计轨迹长度。fillna(0)保证没有轨迹经过的网格不会被统计漏掉。这一步计算量取决于轨迹总顶点数。广州全年轨迹shp如果超过5万条线段可以按月份分组逐月计算再合并否则内存和CPU都会吃紧。网格粒度没有绝对标准建议按分析目的选择。网格大小效果适用场景100米坡度信息保留最好但碎片多单条路线详细查勘200米山体走向清晰热点突出城市徒步活动评估500米平滑但不区分具体小路跨区域绿道对比4.2 从网格热度里提取主要徒步通道有了grid_density后可以用排序快速找到最热门的活动区域。hotspots grid_density.sort_values(seg_len, ascendingFalse) print(hotspots.head(10))只看前10个格子往往发现它们连成一条带状区域这说明白云山或火炉山的主山脊线被频繁踩踏。要提取通道可以把网格按连通性聚合用dissolve把相邻热门格子合并成面再计算面要素与原始轨迹的交叉数量。这一步我一般会把seg_len大于1000米的网格定义为“主力徒步通道”因为200米格子里累计1000米长度意味着至少5次重复通过该区域。4.3 对轨迹点做核密度热力图渔网给出的是离散格子强度视觉上像马赛克。想生成连续热力面可以先把轨迹按固定间距重采样成点再用核密度估算点密度。import numpy as np from scipy.stats import gaussian_kde def resample_line(geom, step10): dist np.arange(0, geom.length, step) pts [geom.interpolate(d) for d in dist] return np.array([(p.x, p.y) for p in pts]) samples np.vstack([ resample_line(g, 10) for g in trails_clean.geometry ]) kde gaussian_kde(samples.T, bw_method0.05)step10代表每10米采一个点避免弯道处点密、直道处点疏的问题。bw_method控制平滑程度数值越大热力面越平滑0.02到0.05适合广州这种中等范围投影坐标。设置成0.05后白云山和火炉山之间不会完全糊成一个整体但也不会出现单条线独占热力尖峰。得到kde对象后在网格范围内生成坐标矩阵再调用kde(points)就能得到每个像素的密度值。这个方法不依赖ArcGIS的密度分析工具适合写进自动化流水线。5. 导出txt、KML并批量出图交付清洗和统计做完了最终交付对象可能是不懂GIS的同事也可能是微信小程序前端。轨迹shp不能直接在网页或导航软件里打开需要导出成txt坐标序列和KML格式同时批量输出按月份拆分的示意图。5.1 shp转TXT输出经度、纬度坐标序列部分户外小程序只接受带经纬度的文本文件。导出前先确认当前几何处于EPSG:4326否则输出的是米制坐标。trails_wgs trails_clean.to_crs(EPSG:4326) frames [] for idx, row in trails_wgs.iterrows(): if row.geometry.geom_type ! LineString: continue pts list(row.geometry.coords) tmp pd.DataFrame({lon: [p[0] for p in pts], lat: [p[1] for p in pts]}) tmp[name] row.get(name, ftrail_{idx}) frames.append(tmp) out_txt pd.concat(frames, ignore_indexTrue) out_txt.to_csv(trails_2020_points.txt, indexFalse)输出的是每一行一个坐标点。如果后续要按线路拆分可以用name字段作为过滤键。geom_type不是LineString时直接跳过避免MultiLineString的坐标遍历炸掉脚本。5.2 用GDAL把轨迹shp转KMLKML的坐标顺序是经度、纬度而且要求使用WGS84坐标系。最稳妥的转换方式是直接用ogr2ogr命令。ogr2ogr -f KML -t_srs EPSG:4326 \ -lco COORDINATE_PRECISION7 \ trails_2020.kml 广州市户外运动轨迹-2020年徒步轨迹.shp-t_srs EPSG:4326在输出阶段做重投影保证所有轨迹都转成经纬度。COORDINATE_PRECISION7允许保留小数点后7位坐标精度约1厘米对轨迹展示足够。转换完成后用QGIS或Google Earth重新打开确认没有跨日期线断线问题。5.3 按月份批量出图并验证输出批量出图用Matplotlib最直接先把所有轨迹画成浅灰背景再高亮当前月份。import matplotlib.pyplot as plt for month, group in trails_clean.groupby(month): fig, ax plt.subplots(figsize(10, 10)) trails_clean.geometry.plot(axax, colorlightgray, linewidth0.5) group.geometry.plot(axax, color#e34a33, linewidth0.8) ax.set_title(fGuangzhou Hiking Trails 2020-{month:02d}) ax.set_axis_off() fig.savefig(fhike_2020_{month:02d}.png, dpi150) plt.close(fig)浅灰背景是所有轨迹红色高亮是当月轨迹这样月与月之间的活动范围一眼就能看出变化。设置ax.set_axis_off()是为了去掉坐标轴刻度避免出图时带上不必要的地理网格线。出完图后最后做一个坐标范围验证。用GeoPandas重新读回KML检查边界是否还在广州范围内。check gpd.read_file(trails_2020.kml) print(check.crs) print(check.total_bounds)total_bounds前两位约113.7和114.5后两位约22.5和23.5说明坐标转换和导出没有发生偏移。如果出现负值或大数多半是原始shp坐标系声明错误需要回到第2章重新检查CRS而不是修正导出的KML。本文还有配套的精品资源点击获取 轨迹SHP数据处理全攻略：从坐标清洗到热度分析
简介面向户外运动研究、城市规划与智慧城市等领域这份GIS数据集提供了广州市2020年徒步活动的完整轨迹记录。数据以标准Shapefile格式存储包含轨迹ID、徒步距离、运动速度、采集日期等核心属性可广泛应用于市民徒步行为分析、热门路线识别、运动设施需求评估等场景适合地理信息数据分析师、城市规划师及运动科学方向的研究者使用。压缩包共7个文件主文件为shp几何数据配套的dbf属性表保存了全部轨迹字段prj文件定义了坐标系统xml元数据提供了数据集说明sbn、sbx、shx索引文件则用于加速空间查询与读取包体约649.82MB。目前已有265人学习或下载。每个轨迹ID可区分不同徒步活动距离与速度属性为运动强度分析提供依据直接使用可省去自行采集GPS轨迹的繁琐流程快速开展徒步空间分布、区域运动资源配置等研究结合天气、人口密度等多源数据还能进行多因素深度挖掘为城市公共空间优化与公共健康决策提供数据支撑。1. 把一份轨迹shp当资产而不是图层收到“广州市户外运动轨迹-2020年徒步轨迹shp”这类数据时第一反应通常是拖进ArcMap看一眼但几秒钟后就会发现满屏都是密密麻麻的折线属性表里全是拼音缩写根本分不清哪条是白云山、哪条是火炉山。这份数据的价值也不在“能显示出来”而在回答“哪段路被走得最多”“哪些线路一年里重复率最高”“热度分布围绕什么山体展开”这类实际问题。shp本身不是轨迹专用格式它把带有时间语义的轨迹压成了二维线要素属性字段还经常混着起点终点所以直接分析会连续踩坐标系、拓扑和字段语义三个坑。下面这套处理流程专为轨迹型shp设计用GeoPandas做读取、清洗和统计用GDAL和QGIS做导出与批量出图不依赖ArcGIS也不需要企业级GIS许可适合户外赛事组织者、城市绿道评估和路网优化方向的工程师。2. 先建立坐标系与字段映射别急着画线轨迹shp最常见的问题不是几何错误而是坐标系声明缺失或字段单位不一致。打开要素前先问自己三个问题坐标是经纬度还是投影米每条线是一条完整的徒步记录还是多个线段叠加属性表里的时间字段能不能直接解析这三件事不解决后面的抽稀、测长、渔网统计全部会失真。2.1 用GeoPandas读取并核对轨迹shp的几何类型先用一段最小化代码把shp的关键信息打印出来这一步能避免后面80%的错误。import geopandas as gpd trails gpd.read_file(广州市户外运动轨迹-2020年徒步轨迹.shp) print(trails.crs) print(trails.geom_type.unique()) print(trails.columns.tolist()) print(trails.head(3))geom_type.unique()会给出LineString或MultiLineString两种典型结果。如果是后者后续length属性和坐标遍历都需要按多部件逐段处理。crs输出None意味着shp缺少.prj文件GeoPandas会按无投影坐标读入此时算出的长度是十进制度数而不是米。columns.tolist()用来确认字段名里是否有start_time、end_time、name、length如果没有时间字段后面的时序拆分就没法做。拿到字段名后还要看一眼坐标范围判断数据是否真的落在广州附近。print(trails.total_bounds)广州的经纬度范围大致是东经113.0至114.5北纬22.5至24.0。如果total_bounds的小数点前是六位数比如[380000, 2600000, 390000, 2620000]说明这份shp已经做过投影不能直接按WGS84处理。2.2 统一坐标系先声明EPSG:4326再投到米制轨迹采集设备输出的原生坐标通常是WGS84经纬度ESRI的shp组件会把投影信息写进.prj但单独压缩传递时经常漏掉这个文件。遇到crs为None的情况不要盲目使用to_crs而是先声明。if trails.crs is None: trails trails.set_crs(EPSG:4326) else: print(已有坐标系, trails.crs)set_crs和to_crs的区别一定要分清前者只是给坐标贴上坐标系标签不改变数值后者才是真正做投影变换。如果坐标范围显示是投影坐标我却把它声明成4326后面的所有长度都会错得离谱。稳妥的判断方式是看total_bounds的绝对量级经纬度在0到180之间米制在10万到上千万之间。测距和做渔网统计前需要转换到适合广州的投影坐标系。trails_m trails.to_crs(EPSG:32649) trails_m[length_m] trails_m.geometry.length print(trails_m[length_m].describe())EPSG:32649是WGS 84 / UTM zone 49N广州恰好落在中央经线117°W附近长度变形很小适合做距离计算和栅格统计。不同坐标系各有用处不是所有米制坐标都适合算长度。EPSG名称单位主要用途4326WGS 84度原始轨迹存储、KML/TXT交换、GPS设备输出32649WGS 84 / UTM zone 49N米长度测量、缓冲区、渔网统计、核密度分析3857Web墨卡托米在线底图叠加显示不适合面积和距离计算坐标转换后再看length_m的统计值如果发现中位数只有几米说明原始shp里的线已经被打断需要先按轨迹ID合并而不是做分析。2.3 时间字段解析把shp变成可按时序分析的数据轨迹shp属性表里常见时间字段有start_time、end_time、date。这些字段在Excel里看着正常读进GeoPandas后可能还是字符串分词需要显式转换。import pandas as pd trails[start_dt] pd.to_datetime( trails[start_time], errorscoerce ) trails[year] trails[start_dt].dt.year trails[month] trails[start_dt].dt.month print(trails[month].value_counts().sort_index())errorscoerce表示解析失败时置为NaT不会因为个别脏数据中断整个脚本。解析完成后马上统计month分布可以判断这条shp是否真的覆盖2020年全年还是只集中在下半年。如果字段原本是时间戳数字比如1609430400需要先除以86400再转日期不能直接丢给pd.to_datetime。时间字段解析完成才谈得上清理几何。很多轨迹shp在采集时没有做航点抽稀一条白云山环线可能包含上万个顶点几何修复的优先级反而高于抽稀。3. 清洗轨迹shp几何修复、抽稀和字段拆分轨迹shp的几何脏点通常有三类自相交导致is_validFalse大量冗余顶点拖慢绘制和计算同一条轨迹被GPS中断拆成多个要素。这三类问题不解决后续聚合统计的结果都会偏高或偏低。清洗的目标不是把几何改得多整齐而是让它满足后续“按线聚合、按米求和”的基本前提。3.1 修复无效几何与自相交线段用Shapely自带的make_valid可以一次性处理大部分几何问题。from shapely.validation import make_valid from shapely.ops import linemerge trails[geometry] trails.geometry.apply(make_valid) def to_single_line(geom): if geom.geom_type MultiLineString: merged linemerge(geom) return merged if merged.geom_type LineString else geom return geom trails[geometry] trails.geometry.apply(to_single_line)make_valid会把自相交的线拆成多个部件所以调用后要用linemerge把首尾相接的部件重新合并成一条完整线。linemerge只处理端点完全重合的情况如果GPS断点之间有几米空隙接口会保持MultiLineString这一步也算变相检测出轨迹中断。修复完成后检查无效几何比例print(无效比例, round((~trails.is_valid).mean(), 4))如果无效比例超过1%说明数据源质量较差后面统计网格时要把无效要素先过滤掉否则overlay可能报错。3.2 抽稀用simplify控制坐标点密度轨迹shp的顶点往往每隔一两米就有一个出图和计算都容易被数据量拖垮。道格拉斯-普克算法的目的不是简化轮廓而是去掉对几何形状影响小于阈值的点。trails_m[geom_simp] trails_m.geometry.simplify( tolerance2.0, preserve_topologyTrue ) trails_clean trails_m.set_geometry(geom_simp).drop(columns[geometry])tolerance单位与坐标系一致。前面我们用了EPSG:32649所以2.0代表2米。容差越大顶点越少但也越容易把急转弯拉直。tolerance适用场景0.5米保留栈道拐角、观景台停留点1.0米普通山地徒步轨迹2.0米广州山体绿道、城市步道5.0米公里级热度底图不适合作精确路线抽稀后记得重新计算长度因为简化会稍微缩短线长。若原始轨迹点非常密集抽稀前后总长度差异应小于1%差异过大说明容差设得太大。3.3 缺失字段拆分与分组输出清洗完几何后很多场景需要按月份或路线名称输出子集。常见的做法是按month字段分组逐组写出独立的shp文件。for month, group in trails_clean.groupby(month): group.to_file( ftrails_2020_{month:02d}.shp, encodingutf-8 ) print(month, len(group), group[length_m].sum())这里的encodingutf-8只影响.dbf属性表。如果后续用ArcGIS打开老版本读UTF-8中文可能乱码一般建议在ArcGIS里也设置代码页。实际上更稳的交付方案是保留字段名用拼音或英文把中文名称放到name字段值里。字段拆分还可以按轨迹距离过滤。比如只保留长度大于1公里的轨迹用来过滤掉GPS漂移产生的碎线。trails_clean trails_clean[trails_clean[length_m] 1000]轨迹清洗阶段不需要过多纠结“线是否光滑”重点是让每个要素都能代表一次独立徒步。下一步要做的是把这些线变成可以排序的密度指标。4. 用渔网和核密度估算徒步热点路线是否热门不能只靠“看起来线多”来判断。把轨迹shp切割到规则网格里统计每个网格内轨迹线的累计长度比目测叠加更客观。这一节先用渔网做分段统计再用核密度做连续热力面。4.1 在投影坐标系下建立200米渔网渔网大小直接影响统计粒度。广州山体徒步轨迹的路线宽度通常在1到3米网格设到100米会切碎同一条山脊线设到500米又分不清白云山和火炉山的边界。200米是折中值。from shapely.geometry import box xmin, ymin, xmax, ymax trails_clean.total_bounds cellsize 200 cols int((xmax - xmin) // cellsize) 1 rows int((ymax - ymin) // cellsize) 1 grid_polys [] for i in range(cols): for j in range(rows): x0 xmin i * cellsize y0 ymin j * cellsize grid_polys.append(box(x0, y0, x0 cellsize, y0 cellsize)) grid gpd.GeoDataFrame(geometrygrid_polys, crstrails_clean.crs) grid[grid_id] range(len(grid))这段代码根据轨迹总范围生成一个全覆盖规则网格每个格子都是独立的矩形要素。grid_id是后续聚合的关键索引不要省略。接下来做线和网格的相交计算。inter gpd.overlay(grid, trails_clean[[geometry]], howintersection) inter[seg_len] inter.geometry.length density ( inter.groupby(grid_id)[seg_len] .sum() .reset_index() ) grid_density grid.merge(density, ongrid_id, howleft) grid_density[seg_len] grid_density[seg_len].fillna(0)overlay会把落在多个网格里的轨迹线逐一分割seg_len是每段在网格内的实际长度。按grid_id汇总后grid_density就带上了每个网格的累计轨迹长度。fillna(0)保证没有轨迹经过的网格不会被统计漏掉。这一步计算量取决于轨迹总顶点数。广州全年轨迹shp如果超过5万条线段可以按月份分组逐月计算再合并否则内存和CPU都会吃紧。网格粒度没有绝对标准建议按分析目的选择。网格大小效果适用场景100米坡度信息保留最好但碎片多单条路线详细查勘200米山体走向清晰热点突出城市徒步活动评估500米平滑但不区分具体小路跨区域绿道对比4.2 从网格热度里提取主要徒步通道有了grid_density后可以用排序快速找到最热门的活动区域。hotspots grid_density.sort_values(seg_len, ascendingFalse) print(hotspots.head(10))只看前10个格子往往发现它们连成一条带状区域这说明白云山或火炉山的主山脊线被频繁踩踏。要提取通道可以把网格按连通性聚合用dissolve把相邻热门格子合并成面再计算面要素与原始轨迹的交叉数量。这一步我一般会把seg_len大于1000米的网格定义为“主力徒步通道”因为200米格子里累计1000米长度意味着至少5次重复通过该区域。4.3 对轨迹点做核密度热力图渔网给出的是离散格子强度视觉上像马赛克。想生成连续热力面可以先把轨迹按固定间距重采样成点再用核密度估算点密度。import numpy as np from scipy.stats import gaussian_kde def resample_line(geom, step10): dist np.arange(0, geom.length, step) pts [geom.interpolate(d) for d in dist] return np.array([(p.x, p.y) for p in pts]) samples np.vstack([ resample_line(g, 10) for g in trails_clean.geometry ]) kde gaussian_kde(samples.T, bw_method0.05)step10代表每10米采一个点避免弯道处点密、直道处点疏的问题。bw_method控制平滑程度数值越大热力面越平滑0.02到0.05适合广州这种中等范围投影坐标。设置成0.05后白云山和火炉山之间不会完全糊成一个整体但也不会出现单条线独占热力尖峰。得到kde对象后在网格范围内生成坐标矩阵再调用kde(points)就能得到每个像素的密度值。这个方法不依赖ArcGIS的密度分析工具适合写进自动化流水线。5. 导出txt、KML并批量出图交付清洗和统计做完了最终交付对象可能是不懂GIS的同事也可能是微信小程序前端。轨迹shp不能直接在网页或导航软件里打开需要导出成txt坐标序列和KML格式同时批量输出按月份拆分的示意图。5.1 shp转TXT输出经度、纬度坐标序列部分户外小程序只接受带经纬度的文本文件。导出前先确认当前几何处于EPSG:4326否则输出的是米制坐标。trails_wgs trails_clean.to_crs(EPSG:4326) frames [] for idx, row in trails_wgs.iterrows(): if row.geometry.geom_type ! LineString: continue pts list(row.geometry.coords) tmp pd.DataFrame({lon: [p[0] for p in pts], lat: [p[1] for p in pts]}) tmp[name] row.get(name, ftrail_{idx}) frames.append(tmp) out_txt pd.concat(frames, ignore_indexTrue) out_txt.to_csv(trails_2020_points.txt, indexFalse)输出的是每一行一个坐标点。如果后续要按线路拆分可以用name字段作为过滤键。geom_type不是LineString时直接跳过避免MultiLineString的坐标遍历炸掉脚本。5.2 用GDAL把轨迹shp转KMLKML的坐标顺序是经度、纬度而且要求使用WGS84坐标系。最稳妥的转换方式是直接用ogr2ogr命令。ogr2ogr -f KML -t_srs EPSG:4326 \ -lco COORDINATE_PRECISION7 \ trails_2020.kml 广州市户外运动轨迹-2020年徒步轨迹.shp-t_srs EPSG:4326在输出阶段做重投影保证所有轨迹都转成经纬度。COORDINATE_PRECISION7允许保留小数点后7位坐标精度约1厘米对轨迹展示足够。转换完成后用QGIS或Google Earth重新打开确认没有跨日期线断线问题。5.3 按月份批量出图并验证输出批量出图用Matplotlib最直接先把所有轨迹画成浅灰背景再高亮当前月份。import matplotlib.pyplot as plt for month, group in trails_clean.groupby(month): fig, ax plt.subplots(figsize(10, 10)) trails_clean.geometry.plot(axax, colorlightgray, linewidth0.5) group.geometry.plot(axax, color#e34a33, linewidth0.8) ax.set_title(fGuangzhou Hiking Trails 2020-{month:02d}) ax.set_axis_off() fig.savefig(fhike_2020_{month:02d}.png, dpi150) plt.close(fig)浅灰背景是所有轨迹红色高亮是当月轨迹这样月与月之间的活动范围一眼就能看出变化。设置ax.set_axis_off()是为了去掉坐标轴刻度避免出图时带上不必要的地理网格线。出完图后最后做一个坐标范围验证。用GeoPandas重新读回KML检查边界是否还在广州范围内。check gpd.read_file(trails_2020.kml) print(check.crs) print(check.total_bounds)total_bounds前两位约113.7和114.5后两位约22.5和23.5说明坐标转换和导出没有发生偏移。如果出现负值或大数多半是原始shp坐标系声明错误需要回到第2章重新检查CRS而不是修正导出的KML。本文还有配套的精品资源点击获取