简介这份1950—2018年中国台风Shapefile数据集面向气象、气候、GIS及灾害风险管理领域的研究者与学习者提供中国沿海台风历史记录的空间化整理可用于台风频率、强度、登陆地点分布、路径特征和时空关联等分析。压缩包共10个文件约1005KB包含shp主数据、shx索引、dbf属性表、prj投影及xml元数据等配套类型构成可直接在ArcGIS、QGIS中打开的点、线两图层点图层对应台风登陆位置线图层记录移动轨迹数据覆盖半个多世纪。已有1785人学习下载说明其具有一定参考价值。读者可据此复现历年台风路径可视化统计年度频次与强度变化识别高风险岸段并结合海平面等资料开展风暴潮影响模拟为城市规划与防灾减灾研究提供历史依据。1. 拿到《1950-2018中国台风shp数据.rar》先急着在GIS里打开看看往往看不出结论手头有一份 1950-2018 年的中国台风 shp 数据多数人的第一反应是解压后用 QGIS 双击打开结果几万个点密密麻麻地铺在海面上放大缩小半天也看不出规律。问题并不在数据本身台风路径本质上是「带时间戳的点序列」shapefile 只负责存点时间、强度、气压这些信息都记录在属性表里不把属性表和几何信息一起看就只能看到一张没有意义的点云图。这份资源的实际价值是把 68 年的台风路径点统一成一套可被 GIS 和 Python 直接消费的矢量结构属性字段规整、坐标系明确省掉了从几十个年度文本文件里拼字段、对表头、补缺失值的体力活。它适合三类人做风电场、光伏、输电线路选址需要评估极端风况的工程人员做沿海城市内涝、风暴潮、保险理赔建模的分析岗以及要复现台风路径可视化、给前端提供轨迹数据的开发者。如果你只是想要一张能交差的图用公开的在线路径图就足够了但如果你要做空间叠加、落区统计、缓冲区分析这套数据是能直接放到分析流程中使用的。2. 解压之后的文件组shp 只是一组文件中的入口文件2.1 .shp、.shx、.dbf、.prj、.cpg 各管什么很多人在解压这个压缩包时是随手操作的把 .shp 单独拷到一个文件夹里剩下几个同名的兄弟文件留在原地等打开时弹出一个「缺少 dbf 或 shx」的报错才发现问题出在解压这一步。shapefile 虽然名字里带 file但它其实是一组文件共同构成的格式.shp 存几何点、线、面坐标.shx 是几何索引.dbf 是属性表本质是个 dBase 表.prj 存坐标系定义.cpg 声明 .dbf 的字符编码.sbn/.sbx 是空间索引可有可无删了不影响读取只是空间查询会慢。所以解压后要做的第一件事是新建一个纯英文路径的工作目录把同名文件整组拷过去中文路径和空格在 fiona、GDAL 这类底层库上偶尔会遇到玄学问题血泪经验是能用英文路径就别挑战。如果压缩包里是按年份或按台风拆了多个 shp先用目录结构看一遍再决定是合并还是抽样。# 建工作目录避免中文路径和空格 mkdir -p ~/work/typhoon/data cd ~/work/typhoon # 解压后看一眼文件组是否完整 ls -lh data/ # 期望看到同名的一组.shp .shx .dbf .prj .cpg这段命令没有技术含量但它把后面 80% 的「读不出来」问题挡在了门外。重点看两件事有没有 .prj决定坐标系能不能自动识别有没有 .cpg决定属性表中文会不会变问号。两个都没有也可以继续但读取时要手工指定第三章会写。2.2 属性表字段用一行代码看清全部字段不同来源的台风路径数据字段名可能是中文、拼音或英文缩写绝不能照着别人的脚本硬套。常见做法是先打印列名和前几行再决定字段映射。下面这张表是这类数据里最典型的几类字段实际名字以你打开看到的为准典型字段类型含义在分析里的用法编号文本或长整型台风年度编号通常「年份 两位序号」分组、去重计数、跨表关联时间文本或日期观测时刻多数是 6 小时一个点排序、断点检测、趋势统计经度 / 纬度双精度地理坐标下的度值已由 geometry 承载属性里常冗余一份中心气压整型单位 hPa强度辅助判断、极值筛选近中心最大风速双精度单位 m/s强度分级的主字段强度等级文本热带低压 / 热带风暴 / 台风 / 强台风 / 超强台风直接筛选、配色分类序号整型该台风内部的点序时间字段缺失时的排序兜底字段看明白后先想清楚你的问题落在哪一列问「这条路径经过哪些地方」看几何问「有多强」看风速和气压问「哪年最多」看编号配合时间。「序号」这一列千万别忽略有些数据的「时间」只有日期没有时刻同一天有好几个点这时候只能靠序号还原顺序否则轨迹会来回折返。提示如果字段名是英文缩写且看不懂不要猜去看随包附带的说明文件没有说明文件就按数值范围倒推风速通常在 1070 之间气压在 8801010 之间经纬度在 100180 和 050 之间。2.3 坐标系判定先看 .prj再看 total_bounds坐标系判断错后面全是无用功。判断顺序只有两步第一步看有没有 .prj有的话用文本编辑器打开是GEOGCS[GCS_WGS_1984...]这类就是地理坐标系EPSG:4326 或等价的 CGCS2000第二步看数值范围total_bounds返回的 x 在 100180、y 在 050 之间说明经纬度没被投影过直接当 4326 用就行如果 x 是六位数、y 是七位数说明已经投影过多半是某个高斯-克吕格或 UTM 分带这时候必须找到带号找不到就别硬算距离。一个容易被忽略的点经纬度顺序颠倒。有些数据在导出时把 x/y 写反了几何上表现为所有点挤在南极附近或者跑到高纬度total_bounds一眼就能发现y 出现负数或 x 出现 050 的数。这种错不是坐标系错是几何本身错了读进来做一次坐标交换即可。3. 用 Python 把路径点读进来GeoPandas 与 pyshp 两条读法3.1 GeoPandas五行代码看清数据全貌GeoPandas 是这类数据最省事的选择底层走 fiona/pyogrio.prj 和 .cpg 会自动读取中文属性基本不用额外处理。import geopandas as gpd # encoding 只在 .cpg 缺失时才需要手工指定 gdf gpd.read_file(data/typhoon_1950_2018.shp, encodingutf-8) print(gdf.crs) # 看坐标系None 说明没有 .prj print(len(gdf)) # 点总数量级通常在几万 print(gdf.columns.tolist()) # 字段清单对照 2.2 的表格 print(gdf.total_bounds) # [minx, miny, maxx, maxy]判断是否投影过 print(gdf.head(3).drop(columnsgeometry))total_bounds是这套流程里最该养成的习惯它在你写任何分析代码之前就告诉你数据的空间范围对不对、坐标系是不是预期。如果 minx 是 100 上下、maxy 是 50 上下说明是经纬度如果 maxy 出现 3000000 这种数就是投影坐标后面所有涉及距离的计算都得统一到同一套投影下再比。读进来之后建议立刻做一次体检几何是否为空、是否有效、编号是否有空值。这三项里任何一项出问题后面的分组和空间连接都会静默丢数据。gdf gdf[gdf.geometry.notna() gdf.geometry.is_valid] gdf gdf.dropna(subset[编号]) print(gdf[编号].nunique()) # 台风个数和年份跨度对得上才算合理3.2 pyshp不装 GDAL 也能读适合流式处理有些机器上装 GDAL 会卡在编译环节pyshp 是纯 Python 实现pip install pyshp一把过读取方式更底层适合只要属性表、或者要逐条流式处理避免一次性载入内存的场景。import shapefile # 没有 .cpg 时手工指定编码中文数据常见 gbk / gb18030 reader shapefile.Reader(data/typhoon_1950_2018.shp, encodinggbk) fields [f[0] for f in reader.fields[1:]] # fields[0] 是 DeletionFlag要跳过 print(fields) print(reader.bbox) # 同 total_bounds print(reader.shapeTypeName) # POINT / POLYLINE决定后续处理方式 # 逐条取适合大数据量下的流式统计 for sr in reader.iterShapeRecords(): rec dict(zip(fields, sr.record)) lon, lat sr.shape.points[0] # 点要素只有一组坐标 print(rec[编号], rec[时间], lon, lat) breakshapeTypeName这一步不能省如果数据是POINT说明一行一个时刻路径要靠编号时间重组如果已经是POLYLINE那每个要素就是一整条路径属性表里存的往往是路径整体的统计值最低气压、最大风速处理方式完全不同。这两种情况在「台风 shp」这个关键词下都能搜到拿到手先确认别默认。3.3 从散点到轨迹线按编号排序并在缺测处打断散点做不了长度、方向、加速度这类分析必须重组成线。核心逻辑就三步按编号分组、组内按时间排序、相邻点连成 LineString。这里有个必须处理的边界原始数据里存在缺测两个点之间可能隔了 12 小时甚至更久直接连起来会画出一条笔直穿越陆地的假路径把它当作真实路径去算落区结果就全错了。import pandas as pd import geopandas as gpd from shapely.geometry import LineString gdf gpd.read_file(data/typhoon_1950_2018.shp, encodingutf-8) gdf gdf[gdf.geometry.notna() gdf.geometry.is_valid].copy() # 时间字段名按实际替换解析失败置 NaT后面直接丢掉 gdf[dt] pd.to_datetime(gdf[时间], errorscoerce) gdf gdf.dropna(subset[dt]).sort_values([编号, dt]) recs [] for tid, g in gdf.groupby(编号): g g.sort_values(dt) if len(g) 2: continue # 相邻点时间差超过 6 小时视为缺测重新分段 gap_h g[dt].diff().dt.total_seconds().div(3600).fillna(0) g g.assign(seg(gap_h 6).cumsum()) for seg_id, sg in g.groupby(seg): if len(sg) 2: continue recs.append({ 编号: tid, 段号: int(seg_id), 起点: sg[dt].iloc[0], 终点: sg[dt].iloc[-1], 点数: len(sg), geometry: LineString(sg.geometry.tolist()), }) tracks gpd.GeoDataFrame(recs, geometrygeometry, crsgdf.crs) tracks.to_file(out/tracks.gpkg, driverGPKG)参数上只有两个地方需要你根据实际情况调整6 小时这个阈值如果原始数据是 3 小时一次就改成 3如果发现断点特别多先回头检查时间解析是不是有一半失败了。crsgdf.crs必须带上否则写出去的文件没有坐标系下次打开又得重判。落盘用 GPKG 而不是 shp是因为 GPKG 支持长字段名、UTF-8 属性和一次存多图层中间结果用它更省事最终交付再另说。4. 实际分析的三类场景落区统计、影响缓冲与强度筛选4.1 与省级边界做空间连接统计「到过哪些省」最常见的需求是「这条路径影响了哪些省级行政区、各有多少次」。做法是把点与行政区面做空间连接within谓词只保留落在面内的点。要注意这里统计的是「路径点落在陆地上的次数」台风大部分时间在海上能匹配上的点本身就少所以这个结果适合叫「登陆/内陆影响统计」不能当成影响范围。import geopandas as gpd pts gpd.read_file(data/typhoon_1950_2018.shp, encodingutf-8).to_crs(4326) prov gpd.read_file(data/admin_province.shp).to_crs(4326) # 新版 geopandas 用 predicate老版是 opwithin joined gpd.sjoin(pts, prov[[省名, geometry]], howinner, predicatewithin) # 按台风去重计数同一个台风在同一个省可能有很多个点 stat (joined.groupby(省名)[编号] .nunique() .sort_values(ascendingFalse) .reset_index(name台风次数)) print(stat.head(10))sjoin的两个参数决定结果对不对howinner只留匹配上的predicate选within点完全在面内还是intersects点接触即算在海陆边界附近的点上会有差异。另外别拿统计结果直接当结论先核对一两条已知路径的首末位置和编号确认连接没串行再往下写报告。4.2 影响半径缓冲区不要在 4326 上直接 buffer「影响范围」这类需求通常用缓冲区近似把每个路径点按固定半径比如 200 公里做缓冲再与目标要素做叠加。这里最容易踩的坑是把经纬度坐标直接buffer(200000)在 4326 下缓冲区单位是「度」200000 度的圆能罩住整个地球结果看起来「全都覆盖到了」实际上是一个彻底错误的结论。import geopandas as gpd pts gpd.read_file(data/typhoon_1950_2018.shp, encodingutf-8).to_crs(4326) # 先投影到等积/等距坐标系再 buffer单位变成米 pts_m pts.to_crs(6933) # 全球等积够用且不用分带 buf_m pts_m.copy() buf_m[geometry] pts_m.geometry.buffer(200_000) # 200 km buf buf_m.to_crs(4326) # 按台风合并得到每个台风的整体影响面 swath buf.dissolve(by编号).reset_index() swath.to_file(out/swath.gpkg, driverGPKG)投影选择上跨纬度大范围数据我一般用全球等积投影如 EPSG:6933代价是形状略有变形如果只分析某一纬度带的近海区域分带 UTM 更精确。dissolve会把同一台风的多个圆合并成一个面几何运算量大几万点全做会比较慢常见做法是先按强度或年份筛一遍再做只对关心的子集算影响面。4.3 强度分段与年度频次把属性表变成可画图的表工程上真正要交付的往往不是图而是两张表一张是按强度等级的分布一张是逐年频次。强度分级国内常用的阈值大致是热带低压 10.817.1 m/s、热带风暴 17.224.4、强热带风暴 24.532.6、台风 32.741.4、强台风 41.550.9、超强台风 51.0 以上。用pd.cut一次性切好别写一长串 if。import pandas as pd df pd.read_csv(out/points.csv) # 或直接从 gdf 转 bins [0, 10.8, 17.2, 24.5, 32.7, 41.5, 51.0, 100] labels [热低, 热带风暴, 强热带风暴, 台风, 强台风, 超强台风, 异常值] df[等级] pd.cut(df[风速], binsbins, labelslabels, rightFalse) # 逐年台风个数按编号去重不是按点计数 df[year] pd.to_datetime(df[时间], errorscoerce).dt.year yearly df.groupby(year)[编号].nunique()rightFalse表示左闭右开避免 17.2 这种边界值同时落进两个区间。按点计数会得到一个荒唐的结论点多的年份台风就多必须用nunique按编号去重。如果最后的年度曲线出现明显的台阶状突变先怀疑是数据在某个年代更换了观测频次比如从 6 小时加密到 3 小时不是气候突变。5. 排查清单读不出、对不上、算不对的五种常见情况5.1 现象属性表里的中文全变成问号或乱码原因几乎只有两个.dbf 的实际编码不是 UTF-8这类数据很多是 GBK/GB18030 导出的或者 .cpg 文件缺失导致读取端按默认编码猜。解决方式是按顺序试先看有没有 .cpg没有就手工指定encodinggbk仍是乱码换gb18030它的字符集更大能兼容 GBK 但反过来不行两个都不行再用encoding_errorsreplace强行读进来然后人工核对少数几条记录确认字段名没被截断。# 依次尝试哪个不乱码用哪个 for enc in (utf-8, gbk, gb18030): try: g gpd.read_file(data/typhoon_1950_2018.shp, encodingenc) print(enc, g.columns.tolist()[:3]) except Exception as e: print(enc, failed:, e)5.2 现象时间列是一串 4 万多的数字原因是原始表在 Excel 里被存成了日期序列号导出 dbf 时又变成了数值型。Excel 的日期起点是 1899-12-30所以 40000 多对应的就是 2010 年前后。解决方式是用pd.to_datetime带origin和unit参数换算转换后抽样核对几条确认年份和编号前缀能对上。# 数值在 20000~50000 之间基本可以判定是 Excel 序列号 s pd.to_datetime(df[时间], unitD, origin1899-12-30, errorscoerce)如果数字只有 4 位比如 1950 这样的那是纯年份按format%Y解析如果是 14 位长整型按format%Y%m%d%H%M%S解析。三种情况写法完全不同先看数值位数再决定不要一律套同一种。5.3 现象轨迹图上有一条横穿陆地的笔直长线原因就是第 3.3 节说的缺测未打断两个相隔很久的点被直接连成一条线。解决办法是先算时间间隔大于阈值就切段得到的分段轨迹要保留「段号」字段后续做长度统计时按段累加而不是按台风累加。另外有一种更隐蔽的情况同一个编号下混进了重复点时间相同但坐标略有差异这会让线出现小折返去重用drop_duplicates(subset[编号,dt])处理。5.4 现象台风点和底图/行政边界错开几十公里先排除最简单的可能底图是 Web 墨卡托EPSG:3857而数据是 4326两个图层没统一投影就叠在一起。这类偏移在低纬度很小、高纬度很明显。排除了投影不统一之后再怀疑几何本身有问题——在某些导出流程里经纬度被写反了表现是所有点集中在一条不合理的纬度带上。用total_bounds和随手抽三个点对照已知坐标两分钟能定位。5.5 现象几万个点画图卡死或者内存爆几万个点对计算机来说不算什么卡住通常是画法的问题逐点plt.plot加底图、或者对每个点单独做缓冲区。解决思路是聚合——用二维直方图或 hexbin 把点密度先压成网格再画或者先按编号dissolve成线再去画。真要做交互展示导出 GeoJSON 之前用simplify()抽稀容差给 0.010.05 度肉眼几乎看不出差别文件能小一个数量级。6. 进阶技巧把 68 年路径压成一张密度网格并固化成一段可重复运行的脚本路径点图看多了会麻木真正能一眼看出「哪片海域、哪个季度台风最密集」的方式是网格统计。做法是给每个点算一个网格索引统计两种值点数和台风编号去重数。前者反映观测密度受观测频次影响后者反映真实的发生频次报告里用后者更站得住。import numpy as np import pandas as pd import geopandas as gpd gdf gpd.read_file(data/typhoon_1950_2018.shp, encodingutf-8).to_crs(4326) gdf[lon] gdf.geometry.x gdf[lat] gdf.geometry.y step 0.5 # 网格边长度近海分析用 0.5全国尺度用 1.0 gdf[gx] (gdf[lon] // step) * step gdf[gy] (gdf[lat] // step) * step grid (gdf.groupby([gx, gy]) .agg(点数(编号, size), 台风数(编号, nunique)) .reset_index()) # 透视成矩阵直接喂给 imshow 或热力图 mat grid.pivot(indexgy, columnsgx, values台风数).fillna(0)step的选择是个权衡0.5 度在近海能看出路径的走向细节1.0 度适合做全国尺度的展示再细就会因为样本稀疏出现大量空洞看起来像噪声。//取整再乘回去是最省的网格化写法注意它对齐的是 0.5 的整数倍跨 0 度经线时会有半个格子的偏移做全球图时要额外处理。我一般会把这一整套流程固化成三段01_read.py负责读入与体检、统一 CRS、输出中间 GPKG02_build.py负责点转线、断点切分、缓冲区03_stats.py只做统计和出图不碰几何。这样每次数据更新只改第一段后面两段原样跑。中间文件统一用 GPKG 落盘命名带上处理日期出问题能回退到上一步而不是从头再来。阶段输入输出关键校验读取原始 shp 文件组points.gpkgtotal_bounds 合理、编号无空值构建points.gpkgtracks.gpkg / swath.gpkg段数、平均点数与年份跨度匹配统计tracks.gpkg逐年表、落区表编号去重后总数与首段一致这里最容易出现的低级错误是把「点数」当成「台风数」写进结论。有一次我把一份季度统计直接交给对方对方拿去和公开年报对数差了将近十倍回头才发现整条流程里所有的count()都该换成nunique()。从那以后我每次做台风相关的聚合都会强制走一遍「总数核对」先算出数据里有多少个不重复编号再和年份跨度、年均个数做一次量级比对对不上就先停下来查不往下写任何图表。希望帮到你。本文还有配套的精品资源点击获取
简介这份1950—2018年中国台风Shapefile数据集面向气象、气候、GIS及灾害风险管理领域的研究者与学习者提供中国沿海台风历史记录的空间化整理可用于台风频率、强度、登陆地点分布、路径特征和时空关联等分析。压缩包共10个文件约1005KB包含shp主数据、shx索引、dbf属性表、prj投影及xml元数据等配套类型构成可直接在ArcGIS、QGIS中打开的点、线两图层点图层对应台风登陆位置线图层记录移动轨迹数据覆盖半个多世纪。已有1785人学习下载说明其具有一定参考价值。读者可据此复现历年台风路径可视化统计年度频次与强度变化识别高风险岸段并结合海平面等资料开展风暴潮影响模拟为城市规划与防灾减灾研究提供历史依据。1. 拿到《1950-2018中国台风shp数据.rar》先急着在GIS里打开看看往往看不出结论手头有一份 1950-2018 年的中国台风 shp 数据多数人的第一反应是解压后用 QGIS 双击打开结果几万个点密密麻麻地铺在海面上放大缩小半天也看不出规律。问题并不在数据本身台风路径本质上是「带时间戳的点序列」shapefile 只负责存点时间、强度、气压这些信息都记录在属性表里不把属性表和几何信息一起看就只能看到一张没有意义的点云图。这份资源的实际价值是把 68 年的台风路径点统一成一套可被 GIS 和 Python 直接消费的矢量结构属性字段规整、坐标系明确省掉了从几十个年度文本文件里拼字段、对表头、补缺失值的体力活。它适合三类人做风电场、光伏、输电线路选址需要评估极端风况的工程人员做沿海城市内涝、风暴潮、保险理赔建模的分析岗以及要复现台风路径可视化、给前端提供轨迹数据的开发者。如果你只是想要一张能交差的图用公开的在线路径图就足够了但如果你要做空间叠加、落区统计、缓冲区分析这套数据是能直接放到分析流程中使用的。2. 解压之后的文件组shp 只是一组文件中的入口文件2.1 .shp、.shx、.dbf、.prj、.cpg 各管什么很多人在解压这个压缩包时是随手操作的把 .shp 单独拷到一个文件夹里剩下几个同名的兄弟文件留在原地等打开时弹出一个「缺少 dbf 或 shx」的报错才发现问题出在解压这一步。shapefile 虽然名字里带 file但它其实是一组文件共同构成的格式.shp 存几何点、线、面坐标.shx 是几何索引.dbf 是属性表本质是个 dBase 表.prj 存坐标系定义.cpg 声明 .dbf 的字符编码.sbn/.sbx 是空间索引可有可无删了不影响读取只是空间查询会慢。所以解压后要做的第一件事是新建一个纯英文路径的工作目录把同名文件整组拷过去中文路径和空格在 fiona、GDAL 这类底层库上偶尔会遇到玄学问题血泪经验是能用英文路径就别挑战。如果压缩包里是按年份或按台风拆了多个 shp先用目录结构看一遍再决定是合并还是抽样。# 建工作目录避免中文路径和空格 mkdir -p ~/work/typhoon/data cd ~/work/typhoon # 解压后看一眼文件组是否完整 ls -lh data/ # 期望看到同名的一组.shp .shx .dbf .prj .cpg这段命令没有技术含量但它把后面 80% 的「读不出来」问题挡在了门外。重点看两件事有没有 .prj决定坐标系能不能自动识别有没有 .cpg决定属性表中文会不会变问号。两个都没有也可以继续但读取时要手工指定第三章会写。2.2 属性表字段用一行代码看清全部字段不同来源的台风路径数据字段名可能是中文、拼音或英文缩写绝不能照着别人的脚本硬套。常见做法是先打印列名和前几行再决定字段映射。下面这张表是这类数据里最典型的几类字段实际名字以你打开看到的为准典型字段类型含义在分析里的用法编号文本或长整型台风年度编号通常「年份 两位序号」分组、去重计数、跨表关联时间文本或日期观测时刻多数是 6 小时一个点排序、断点检测、趋势统计经度 / 纬度双精度地理坐标下的度值已由 geometry 承载属性里常冗余一份中心气压整型单位 hPa强度辅助判断、极值筛选近中心最大风速双精度单位 m/s强度分级的主字段强度等级文本热带低压 / 热带风暴 / 台风 / 强台风 / 超强台风直接筛选、配色分类序号整型该台风内部的点序时间字段缺失时的排序兜底字段看明白后先想清楚你的问题落在哪一列问「这条路径经过哪些地方」看几何问「有多强」看风速和气压问「哪年最多」看编号配合时间。「序号」这一列千万别忽略有些数据的「时间」只有日期没有时刻同一天有好几个点这时候只能靠序号还原顺序否则轨迹会来回折返。提示如果字段名是英文缩写且看不懂不要猜去看随包附带的说明文件没有说明文件就按数值范围倒推风速通常在 1070 之间气压在 8801010 之间经纬度在 100180 和 050 之间。2.3 坐标系判定先看 .prj再看 total_bounds坐标系判断错后面全是无用功。判断顺序只有两步第一步看有没有 .prj有的话用文本编辑器打开是GEOGCS[GCS_WGS_1984...]这类就是地理坐标系EPSG:4326 或等价的 CGCS2000第二步看数值范围total_bounds返回的 x 在 100180、y 在 050 之间说明经纬度没被投影过直接当 4326 用就行如果 x 是六位数、y 是七位数说明已经投影过多半是某个高斯-克吕格或 UTM 分带这时候必须找到带号找不到就别硬算距离。一个容易被忽略的点经纬度顺序颠倒。有些数据在导出时把 x/y 写反了几何上表现为所有点挤在南极附近或者跑到高纬度total_bounds一眼就能发现y 出现负数或 x 出现 050 的数。这种错不是坐标系错是几何本身错了读进来做一次坐标交换即可。3. 用 Python 把路径点读进来GeoPandas 与 pyshp 两条读法3.1 GeoPandas五行代码看清数据全貌GeoPandas 是这类数据最省事的选择底层走 fiona/pyogrio.prj 和 .cpg 会自动读取中文属性基本不用额外处理。import geopandas as gpd # encoding 只在 .cpg 缺失时才需要手工指定 gdf gpd.read_file(data/typhoon_1950_2018.shp, encodingutf-8) print(gdf.crs) # 看坐标系None 说明没有 .prj print(len(gdf)) # 点总数量级通常在几万 print(gdf.columns.tolist()) # 字段清单对照 2.2 的表格 print(gdf.total_bounds) # [minx, miny, maxx, maxy]判断是否投影过 print(gdf.head(3).drop(columnsgeometry))total_bounds是这套流程里最该养成的习惯它在你写任何分析代码之前就告诉你数据的空间范围对不对、坐标系是不是预期。如果 minx 是 100 上下、maxy 是 50 上下说明是经纬度如果 maxy 出现 3000000 这种数就是投影坐标后面所有涉及距离的计算都得统一到同一套投影下再比。读进来之后建议立刻做一次体检几何是否为空、是否有效、编号是否有空值。这三项里任何一项出问题后面的分组和空间连接都会静默丢数据。gdf gdf[gdf.geometry.notna() gdf.geometry.is_valid] gdf gdf.dropna(subset[编号]) print(gdf[编号].nunique()) # 台风个数和年份跨度对得上才算合理3.2 pyshp不装 GDAL 也能读适合流式处理有些机器上装 GDAL 会卡在编译环节pyshp 是纯 Python 实现pip install pyshp一把过读取方式更底层适合只要属性表、或者要逐条流式处理避免一次性载入内存的场景。import shapefile # 没有 .cpg 时手工指定编码中文数据常见 gbk / gb18030 reader shapefile.Reader(data/typhoon_1950_2018.shp, encodinggbk) fields [f[0] for f in reader.fields[1:]] # fields[0] 是 DeletionFlag要跳过 print(fields) print(reader.bbox) # 同 total_bounds print(reader.shapeTypeName) # POINT / POLYLINE决定后续处理方式 # 逐条取适合大数据量下的流式统计 for sr in reader.iterShapeRecords(): rec dict(zip(fields, sr.record)) lon, lat sr.shape.points[0] # 点要素只有一组坐标 print(rec[编号], rec[时间], lon, lat) breakshapeTypeName这一步不能省如果数据是POINT说明一行一个时刻路径要靠编号时间重组如果已经是POLYLINE那每个要素就是一整条路径属性表里存的往往是路径整体的统计值最低气压、最大风速处理方式完全不同。这两种情况在「台风 shp」这个关键词下都能搜到拿到手先确认别默认。3.3 从散点到轨迹线按编号排序并在缺测处打断散点做不了长度、方向、加速度这类分析必须重组成线。核心逻辑就三步按编号分组、组内按时间排序、相邻点连成 LineString。这里有个必须处理的边界原始数据里存在缺测两个点之间可能隔了 12 小时甚至更久直接连起来会画出一条笔直穿越陆地的假路径把它当作真实路径去算落区结果就全错了。import pandas as pd import geopandas as gpd from shapely.geometry import LineString gdf gpd.read_file(data/typhoon_1950_2018.shp, encodingutf-8) gdf gdf[gdf.geometry.notna() gdf.geometry.is_valid].copy() # 时间字段名按实际替换解析失败置 NaT后面直接丢掉 gdf[dt] pd.to_datetime(gdf[时间], errorscoerce) gdf gdf.dropna(subset[dt]).sort_values([编号, dt]) recs [] for tid, g in gdf.groupby(编号): g g.sort_values(dt) if len(g) 2: continue # 相邻点时间差超过 6 小时视为缺测重新分段 gap_h g[dt].diff().dt.total_seconds().div(3600).fillna(0) g g.assign(seg(gap_h 6).cumsum()) for seg_id, sg in g.groupby(seg): if len(sg) 2: continue recs.append({ 编号: tid, 段号: int(seg_id), 起点: sg[dt].iloc[0], 终点: sg[dt].iloc[-1], 点数: len(sg), geometry: LineString(sg.geometry.tolist()), }) tracks gpd.GeoDataFrame(recs, geometrygeometry, crsgdf.crs) tracks.to_file(out/tracks.gpkg, driverGPKG)参数上只有两个地方需要你根据实际情况调整6 小时这个阈值如果原始数据是 3 小时一次就改成 3如果发现断点特别多先回头检查时间解析是不是有一半失败了。crsgdf.crs必须带上否则写出去的文件没有坐标系下次打开又得重判。落盘用 GPKG 而不是 shp是因为 GPKG 支持长字段名、UTF-8 属性和一次存多图层中间结果用它更省事最终交付再另说。4. 实际分析的三类场景落区统计、影响缓冲与强度筛选4.1 与省级边界做空间连接统计「到过哪些省」最常见的需求是「这条路径影响了哪些省级行政区、各有多少次」。做法是把点与行政区面做空间连接within谓词只保留落在面内的点。要注意这里统计的是「路径点落在陆地上的次数」台风大部分时间在海上能匹配上的点本身就少所以这个结果适合叫「登陆/内陆影响统计」不能当成影响范围。import geopandas as gpd pts gpd.read_file(data/typhoon_1950_2018.shp, encodingutf-8).to_crs(4326) prov gpd.read_file(data/admin_province.shp).to_crs(4326) # 新版 geopandas 用 predicate老版是 opwithin joined gpd.sjoin(pts, prov[[省名, geometry]], howinner, predicatewithin) # 按台风去重计数同一个台风在同一个省可能有很多个点 stat (joined.groupby(省名)[编号] .nunique() .sort_values(ascendingFalse) .reset_index(name台风次数)) print(stat.head(10))sjoin的两个参数决定结果对不对howinner只留匹配上的predicate选within点完全在面内还是intersects点接触即算在海陆边界附近的点上会有差异。另外别拿统计结果直接当结论先核对一两条已知路径的首末位置和编号确认连接没串行再往下写报告。4.2 影响半径缓冲区不要在 4326 上直接 buffer「影响范围」这类需求通常用缓冲区近似把每个路径点按固定半径比如 200 公里做缓冲再与目标要素做叠加。这里最容易踩的坑是把经纬度坐标直接buffer(200000)在 4326 下缓冲区单位是「度」200000 度的圆能罩住整个地球结果看起来「全都覆盖到了」实际上是一个彻底错误的结论。import geopandas as gpd pts gpd.read_file(data/typhoon_1950_2018.shp, encodingutf-8).to_crs(4326) # 先投影到等积/等距坐标系再 buffer单位变成米 pts_m pts.to_crs(6933) # 全球等积够用且不用分带 buf_m pts_m.copy() buf_m[geometry] pts_m.geometry.buffer(200_000) # 200 km buf buf_m.to_crs(4326) # 按台风合并得到每个台风的整体影响面 swath buf.dissolve(by编号).reset_index() swath.to_file(out/swath.gpkg, driverGPKG)投影选择上跨纬度大范围数据我一般用全球等积投影如 EPSG:6933代价是形状略有变形如果只分析某一纬度带的近海区域分带 UTM 更精确。dissolve会把同一台风的多个圆合并成一个面几何运算量大几万点全做会比较慢常见做法是先按强度或年份筛一遍再做只对关心的子集算影响面。4.3 强度分段与年度频次把属性表变成可画图的表工程上真正要交付的往往不是图而是两张表一张是按强度等级的分布一张是逐年频次。强度分级国内常用的阈值大致是热带低压 10.817.1 m/s、热带风暴 17.224.4、强热带风暴 24.532.6、台风 32.741.4、强台风 41.550.9、超强台风 51.0 以上。用pd.cut一次性切好别写一长串 if。import pandas as pd df pd.read_csv(out/points.csv) # 或直接从 gdf 转 bins [0, 10.8, 17.2, 24.5, 32.7, 41.5, 51.0, 100] labels [热低, 热带风暴, 强热带风暴, 台风, 强台风, 超强台风, 异常值] df[等级] pd.cut(df[风速], binsbins, labelslabels, rightFalse) # 逐年台风个数按编号去重不是按点计数 df[year] pd.to_datetime(df[时间], errorscoerce).dt.year yearly df.groupby(year)[编号].nunique()rightFalse表示左闭右开避免 17.2 这种边界值同时落进两个区间。按点计数会得到一个荒唐的结论点多的年份台风就多必须用nunique按编号去重。如果最后的年度曲线出现明显的台阶状突变先怀疑是数据在某个年代更换了观测频次比如从 6 小时加密到 3 小时不是气候突变。5. 排查清单读不出、对不上、算不对的五种常见情况5.1 现象属性表里的中文全变成问号或乱码原因几乎只有两个.dbf 的实际编码不是 UTF-8这类数据很多是 GBK/GB18030 导出的或者 .cpg 文件缺失导致读取端按默认编码猜。解决方式是按顺序试先看有没有 .cpg没有就手工指定encodinggbk仍是乱码换gb18030它的字符集更大能兼容 GBK 但反过来不行两个都不行再用encoding_errorsreplace强行读进来然后人工核对少数几条记录确认字段名没被截断。# 依次尝试哪个不乱码用哪个 for enc in (utf-8, gbk, gb18030): try: g gpd.read_file(data/typhoon_1950_2018.shp, encodingenc) print(enc, g.columns.tolist()[:3]) except Exception as e: print(enc, failed:, e)5.2 现象时间列是一串 4 万多的数字原因是原始表在 Excel 里被存成了日期序列号导出 dbf 时又变成了数值型。Excel 的日期起点是 1899-12-30所以 40000 多对应的就是 2010 年前后。解决方式是用pd.to_datetime带origin和unit参数换算转换后抽样核对几条确认年份和编号前缀能对上。# 数值在 20000~50000 之间基本可以判定是 Excel 序列号 s pd.to_datetime(df[时间], unitD, origin1899-12-30, errorscoerce)如果数字只有 4 位比如 1950 这样的那是纯年份按format%Y解析如果是 14 位长整型按format%Y%m%d%H%M%S解析。三种情况写法完全不同先看数值位数再决定不要一律套同一种。5.3 现象轨迹图上有一条横穿陆地的笔直长线原因就是第 3.3 节说的缺测未打断两个相隔很久的点被直接连成一条线。解决办法是先算时间间隔大于阈值就切段得到的分段轨迹要保留「段号」字段后续做长度统计时按段累加而不是按台风累加。另外有一种更隐蔽的情况同一个编号下混进了重复点时间相同但坐标略有差异这会让线出现小折返去重用drop_duplicates(subset[编号,dt])处理。5.4 现象台风点和底图/行政边界错开几十公里先排除最简单的可能底图是 Web 墨卡托EPSG:3857而数据是 4326两个图层没统一投影就叠在一起。这类偏移在低纬度很小、高纬度很明显。排除了投影不统一之后再怀疑几何本身有问题——在某些导出流程里经纬度被写反了表现是所有点集中在一条不合理的纬度带上。用total_bounds和随手抽三个点对照已知坐标两分钟能定位。5.5 现象几万个点画图卡死或者内存爆几万个点对计算机来说不算什么卡住通常是画法的问题逐点plt.plot加底图、或者对每个点单独做缓冲区。解决思路是聚合——用二维直方图或 hexbin 把点密度先压成网格再画或者先按编号dissolve成线再去画。真要做交互展示导出 GeoJSON 之前用simplify()抽稀容差给 0.010.05 度肉眼几乎看不出差别文件能小一个数量级。6. 进阶技巧把 68 年路径压成一张密度网格并固化成一段可重复运行的脚本路径点图看多了会麻木真正能一眼看出「哪片海域、哪个季度台风最密集」的方式是网格统计。做法是给每个点算一个网格索引统计两种值点数和台风编号去重数。前者反映观测密度受观测频次影响后者反映真实的发生频次报告里用后者更站得住。import numpy as np import pandas as pd import geopandas as gpd gdf gpd.read_file(data/typhoon_1950_2018.shp, encodingutf-8).to_crs(4326) gdf[lon] gdf.geometry.x gdf[lat] gdf.geometry.y step 0.5 # 网格边长度近海分析用 0.5全国尺度用 1.0 gdf[gx] (gdf[lon] // step) * step gdf[gy] (gdf[lat] // step) * step grid (gdf.groupby([gx, gy]) .agg(点数(编号, size), 台风数(编号, nunique)) .reset_index()) # 透视成矩阵直接喂给 imshow 或热力图 mat grid.pivot(indexgy, columnsgx, values台风数).fillna(0)step的选择是个权衡0.5 度在近海能看出路径的走向细节1.0 度适合做全国尺度的展示再细就会因为样本稀疏出现大量空洞看起来像噪声。//取整再乘回去是最省的网格化写法注意它对齐的是 0.5 的整数倍跨 0 度经线时会有半个格子的偏移做全球图时要额外处理。我一般会把这一整套流程固化成三段01_read.py负责读入与体检、统一 CRS、输出中间 GPKG02_build.py负责点转线、断点切分、缓冲区03_stats.py只做统计和出图不碰几何。这样每次数据更新只改第一段后面两段原样跑。中间文件统一用 GPKG 落盘命名带上处理日期出问题能回退到上一步而不是从头再来。阶段输入输出关键校验读取原始 shp 文件组points.gpkgtotal_bounds 合理、编号无空值构建points.gpkgtracks.gpkg / swath.gpkg段数、平均点数与年份跨度匹配统计tracks.gpkg逐年表、落区表编号去重后总数与首段一致这里最容易出现的低级错误是把「点数」当成「台风数」写进结论。有一次我把一份季度统计直接交给对方对方拿去和公开年报对数差了将近十倍回头才发现整条流程里所有的count()都该换成nunique()。从那以后我每次做台风相关的聚合都会强制走一遍「总数核对」先算出数据里有多少个不重复编号再和年份跨度、年均个数做一次量级比对对不上就先停下来查不往下写任何图表。希望帮到你。本文还有配套的精品资源点击获取
知识建站
内容驱动
阅读友好版式
SEO 精细化
线索转化