1950-2018中国台风shp：GeoPandas轨迹重建与避坑指南

简介覆盖1950至2018年中国沿海台风历史记录的Shapefile矢量数据集面向气象与气候研究、灾害风险管理、城市规划及GIS空间分析的从业者与高年级学生用于复盘台风登陆位置与移动路径支撑频率统计、强度分级、路径特征与风暴潮模拟等分析需求。压缩包共10个文件约1005KB内含点图层与线图层两套数据shp存储几何要素dbf保存属性字段shx提供空间索引prj定义坐标参考xml记录元数据导入ArcGIS、QGIS等软件即可直接调用并叠加制图。已有1785人学习下载在相关教学与科研场景中具有一定认可度。近七十年的时间跨度配合点线双图层便于把登陆点与轨迹分段还原读者可据此统计逐年台风数量、比较不同等级台风的路径差异或结合海平面资料开展风险区划为防灾减灾方案与论文写作提供可直接使用的底图素材。1. 从一张点表说起1950-2018中国台风shp到底能干什么很多人拿到一份1950-2018中国台风shp数据第一反应是拖进QGIS看一眼路径图就收工结果真到做登陆频次统计、风场叠加或者灾损关联时才发现时间列是文本、台风编号对不齐、坐标系还是未经处理的地理坐标——图能出分析全错。这套台风shp本质上是一份跨近七十年的台风路径空间底表每个观测时次对应一条台风中心点记录常见字段覆盖编号、名称、时间、经纬度、风速、气压。它适合做沿海台风频次、强度演变、路径聚类的人也适合需要把台风轨迹接进自建GIS分析管线的开发者。它解决的不是出图好不好看而是有没有一份能直接进Python和GIS流程、字段能对齐、时间能排序的台风轨迹底表。下面按结构、读取、重建、避坑、验证五步拆开讲。2. 台风shp的结构解剖文件组、字段与坐标系动手之前先把数据本身看明白比直接写代码重要得多。这一章解决三个问题压缩包里到底有哪几个文件、属性表里可能有哪些字段、坐标系是经纬度还是投影坐标。这三件事任何一件判断错后面轨迹重建和统计分析都会连锁出错。2.1 shapefile文件组别只盯着.shpshapefile 从来不是单个文件而是一组同名文件协同工作。.shp存几何.shx存几何索引.dbf存属性表.prj存坐标系定义.cpg存属性表的字符编码。解压之后第一件事是列目录确认这套文件是否齐全缺.shx会导致索引建不起来缺.dbf就只剩空几何、一个属性字段都读不到。# rar 包用 unrar 解先解到独立目录别污染工作区 unrar x 1950-2018中国台风SHp数据.rar typhoon_data/ # 列出解压结果重点看同名文件组是否齐全 ls -la typhoon_data/unrar x的x表示保留压缩包内的目录结构解到typhoon_data/下便于后续统一路径。列目录时重点看.shp、.shx、.dbf、.prj、.cpg这五类后缀是否成套出现缺.prj说明坐标系信息没带出来后面得手工判定缺.cpg说明编码没标注中文字段大概率要试编码。注意DBF 的字段名最多 10 个字符超长或中文列名会被截断读出来的列名往往和你想象的英文名不一致。别急着在代码里写死列名先探查再映射。2.2 属性表字段解读编号、时间、风压三件套台风路径数据的字段命名没有强制标准不同来源的列名差异很大但语义基本围绕几类核心信息。先用 GeoPandas 把字段和前几行打出来再决定怎么映射比凭经验写死列名稳得多。import geopandas as gpd # 路径换成你解压出来的 .shp 文件名 gdf gpd.read_file(typhoon_data/typhoon.shp, encodingutf-8) print(要素数:, len(gdf)) print(字段列表:, gdf.columns.tolist()) print(几何类型:, gdf.geom_type.unique()) print(gdf.head(3).to_string()) print(坐标系:, gdf.crs)这段代码先回答四个问题一共多少条点记录、列名到底叫什么、几何是点还是线、坐标系有没有带出来。geom_type如果是Point说明这是一张逐时次的中心点表轨迹线要靠后续分组重建如果已经是LineString那说明数据已经做过轨迹合成处理思路要换。head打印时用to_string()避免列被折叠方便逐列核对。常见字段的语义大致如下表实际命名以你探查结果为准不要照抄列名。语义常见命名类型倾向处理要点台风编号ID / TC_ID / 编号字符或数值统一为字符串防前导零丢失台风名称NAME / 名称字符注意中文编码观测时间TIME / DATE / 时间字符或日期先看原样再转 datetime中心经度LON / 经度数值确认是度不是弧度中心纬度LAT / 纬度数值与经度成对校验最大风速WIND / 风速数值确认单位是 m/s 还是 kt中心气压PRES / 气压数值单位一般为 hPa强度等级GRADE / 等级字符可能为空需回填编号、时间、风速气压是后续所有分析的基石。编号决定分组时间决定排序风速气压决定强度分级三者任何一个类型不对轨迹就会乱窜。2.3 坐标系判定地理坐标还是投影坐标坐标系判断错是台风数据处理里最隐蔽的翻车点。地理坐标经纬度的经度范围在 -180 到 180、纬度在 -90 到 90投影坐标的数值往往是几十万到几百万。光看.crs有时还不够因为部分数据.prj缺失、crs为None这时得靠数值范围反推。print(crs:, gdf.crs) # total_bounds 返回 (minx, miny, maxx, maxy)用它判断数值量级 print(范围:, gdf.total_bounds) # 若确认是 WGS84 经纬度需要算距离或面积时再投影到米制 gdf_m gdf.to_crs(EPSG:32650) # UTM 50N按研究区经度带选择 print(投影后范围:, gdf_m.total_bounds)total_bounds的四个值如果落在个位到百位之间基本就是经纬度如果出现几十万、几百万那是投影坐标不能直接当经纬度塞进box()或within()做空间筛选。to_crs(EPSG:32650)把经纬度转到 UTM 50N 带单位变成米适合算路径长度和缓冲区如果只是画图叠加在线底图用EPSG:3857更省事。选哪个投影带取决于研究区经度横跨多个带时老老实实回到经纬度做拓扑判断只在算距离时临时投影。提示台风点如果出现经度几百万、纬度几十万的组合先别怀疑数据错八成是投影坐标。用total_bounds反推比盯着.crs可靠。3. 用GeoPandas读取与筛选从环境到空间子集第二章看清了数据长什么样这一章把它真正读进分析流程并切出你要的研究子集。环境、读取、时间解析、空间筛选四步每一步都有容易忽略的参数。3.1 环境搭建与依赖GeoPandas 依赖 shapely、pyproj、fiona 三个底层库版本不匹配时会出现各种读文件报错最常见的是 shapely 版本和 pyproj 版本打架。用独立 conda 环境隔离别在 base 里反复装卸。conda create -n typhoon python3.10 -y conda activate typhoon pip install geopandas shapely pandas matplotlib pyprojpython3.10是当前和主流 GIS 库兼容性较稳的版本。geopandas会自动带上shapely、pyproj、fiona等依赖但显式再列一遍能避免 pip 解析出过老的版本。装完先跑python -c import geopandas; print(geopandas.__version__)确认能导入再动数据。注意如果读取时报fiona相关的驱动错误优先怀疑fiona和gdal版本冲突重装fiona往往能解别急着改数据。3.2 读取与字段探查读取之后不要立刻进入分析先把几何类型、空几何、时间列原样打印出来。台风数据常见的问题是时间列混了多种格式直接to_datetime会静默产生大量NaT等分组时才发现大半数据消失。import geopandas as gpd import pandas as pd gdf gpd.read_file(typhoon_data/typhoon.shp, encodingutf-8) print(几何类型:, gdf.geom_type.unique()) print(空几何数:, gdf.geometry.isna().sum()) # 先看时间列原样别急着转 print(时间列样例:, gdf[TIME].head(10).tolist()) gdf[dt] pd.to_datetime(gdf[TIME], errorscoerce) print(解析失败行数:, gdf[dt].isna().sum())errorscoerce的作用是把解析不了的值变成NaT而不是抛异常这样你能先统计失败比例。如果失败行数占比很高说明这一列混了多种格式比如一部分是1950/1/1、一部分是1950-01-01 08:00这时要先用正则或字符串切片统一成一种格式再转而不是硬转。空几何数用来确认有没有坐标缺失的点这些点在重建轨迹时必须剔除否则LineString会直接报错。字段映射这一步建议单独写一个字典把探查到的真实列名对到语义名后面代码统一用语义名换一份数据只改字典不改逻辑。col_map { ID: typhoon_id, NAME: typhoon_name, TIME: raw_time, LON: lon, LAT: lat, WIND: wind, PRES: pressure, } gdf gdf.rename(columns{k: v for k, v in col_map.items() if k in gdf.columns})rename前先判断k in gdf.columns避免某个字段在这份数据里不存在时报KeyError。这一步看起来琐碎但它是后面所有代码能稳定跑的前提——不同来源的台风 shp 列名真的可能完全不同。3.3 按年份与空间范围筛选有了可用的时间列和几何就可以切研究子集。按年份切用dt的.dt.year按空间切用shapely.geometry.box构造矩形范围。注意within和intersects的区别within要求点完全落在框内intersects只要相交就算边界点用within会被排除。# 按年份筛先复制避免 SettingWithCopyWarning gdf[year] gdf[dt].dt.year recent gdf[(gdf[year] 2000) (gdf[year] 2018)].copy() # 按空间范围筛box 参数顺序是 (minx, miny, maxx, maxy)单位与 crs 一致 from shapely.geometry import box aoi box(105, 5, 135, 45) # 经纬度范围按研究区改 subset recent[recent.within(aoi)] print(原始:, len(gdf), 按年:, len(recent), 按范围:, len(subset))box(105, 5, 135, 45)的四个参数是经度最小值、纬度最小值、经度最大值、纬度最大值对应研究区的西南角和东北角。如果数据是投影坐标这个box必须换成投影坐标下的范围否则筛选结果为空或错位。打印三步的数量变化是为了确认筛选逻辑符合预期——如果按年之后数量骤降多半是时间解析失败回头查dt的NaT比例。提示筛选后一定.copy()否则后续给子集加列会触发 pandas 的链式赋值警告甚至在旧版本里静默不生效。4. 从散点到路径线台风轨迹重建与强度分段点表本身只能回答“某时刻台风在哪”要回答“台风怎么走”必须按编号分组、按时间排序、再连成线。这一章把散点变成轨迹并按风速做强度分段最后落一张能看的图。4.1 按台风编号分组与时间排序轨迹重建的前提是同一台风的所有点按时间先后排好。排序键必须同时包含编号和时间只按时间排会把不同台风的点交错在一起连出来的线满屏乱窜。# 先剔除时间为空和编号为空的记录它们无法参与分组排序 clean gdf.dropna(subset[dt, typhoon_id]).copy() # 编号统一成字符串去掉可能的前导空格 clean[typhoon_id] clean[typhoon_id].astype(str).str.strip() # 双键排序先按台风编号再按时间 clean clean.sort_values([typhoon_id, dt]).reset_index(dropTrue) print(clean.groupby(typhoon_id).size().describe())groupby(typhoon_id).size().describe()会给出每个台风的观测点数分布最小值如果是 1说明有些台风只有单点记录这类记录连不成线后面要单独处理。astype(str).str.strip()是为了防止编号里混入空格导致同一台风被拆成两组这是很常见的脏数据形态。4.2 生成LineString轨迹分组之后对每组按时间顺序取几何点点数大于等于 2 才能构成线。用shapely.geometry.LineString直接连groupby.apply返回的是以编号为索引的 Series再重新组装成 GeoDataFrame。from shapely.geometry import LineString def to_line(group): pts list(group.geometry) if len(pts) 2: return None return LineString(pts) tracks (clean.groupby(typhoon_id, group_keysFalse) .apply(to_line) .dropna()) track_gdf gpd.GeoDataFrame( {typhoon_id: tracks.index}, geometrytracks.values, crsclean.crs, ) print(轨迹数:, len(track_gdf))group_keysFalse避免apply额外加一层分组键索引导致结构混乱dropna()去掉单点台风返回的None。GeoDataFrame构造时把tracks.index作为编号列tracks.values作为几何crs继承原点表的坐标系——这一步别漏crs否则后续叠加底图会错位。注意groupby.apply在点数很多时较慢几万条以内没问题如果数据量再大改成先按编号取geometry.tolist()再循环构造LineString性能会好一些。4.3 强度分级与颜色映射轨迹线有了接下来给每个点打强度标签。台风强度按近中心最大风速分级阈值用 m/s 是气象上的通行做法但如果数据里的风速单位是节kt必须先乘 0.514 换算否则分级全错。def wind_grade(w): if pd.isna(w): return 未知 if w 17.2: return 热带低压 elif w 24.5: return 热带风暴 elif w 32.7: return 强热带风暴 elif w 41.5: return 台风 elif w 51.0: return 强台风 else: return 超强台风 # 若探查发现风速是 kt先换算clean[wind] clean[wind] * 0.514 clean[grade] clean[wind].apply(wind_grade) print(clean[grade].value_counts())分级阈值对应热带低压、热带风暴、强热带风暴、台风、强台风、超强台风六档是业务上常用的划分。pd.isna(w)兜住空值避免NaN参与比较时全部落到else被误判成超强台风。value_counts()打印各档数量如果某一档数量为 0 或异常多回头查风速单位是否换算过。画图时把轨迹线和强度点叠在一起线用浅灰、点按等级着色一眼就能看出哪段是快速增强、哪段是减弱。import matplotlib.pyplot as plt fig, ax plt.subplots(figsize(10, 8)) track_gdf.plot(axax, colorlightgray, linewidth0.6) clean.plot(axax, columngrade, markersize6, legendTrue, categoricalTrue) ax.set_title(Typhoon Tracks and Intensity) plt.tight_layout() plt.show()categoricalTrue让分级按类别着色而不是连续色带图例更可读。markersize6在近七十年全量数据下偏大如果点太密可以调到 2 到 3或者先按年份筛完再画。5. 台风shp处理避坑五条血泪经验这一章单独列五条实际处理中最常翻车的问题每条按现象、原因、解决三段走。这些问题不会在第一次跑代码时全部暴露往往在数据量放大或换一份来源后才冒出来。5.1 中文属性乱码问号与方块现象读进来的名称字段全是问号、方块或乱码英文编号正常。 原因DBF 的属性表编码和读取时指定的encoding不一致.cpg文件缺失或标注错误时最明显GBK 编码的数据用 UTF-8 读就会乱。 解决先看解压目录有没有.cpg有就按它标的编码读没有就依次试encodinggbk和encodingutf-8用gdf.head()看名称列是否正常。别用latin-1硬读那只会把乱码变成另一种乱码。5.2 台风编号前导零丢失现象编号0123读进来变成123分组后轨迹数量对不上或者同一台风被拆成两条。 原因DBF 里编号是数值型字段GeoPandas 读成整型前导零被吃掉。 解决读取后立刻astype(str)需要补齐位数就str.zfill(4)。更稳的做法是在read_file后用dtype{ID: str}显式指定但字段名要先探查确认。5.3 时间列跨格式导致大量NaT现象pd.to_datetime之后isna().sum()占比很高分组时大量点消失。 原因同一列混了1950/1/1、1950-01-01、1950-01-01 08:00多种写法甚至带时区后缀。 解决先用errorscoerce转一轮把NaT的行单独抽出来看原始字符串长什么样再用字符串替换统一分隔符最后指定format精确解析。不要指望to_datetime自动兜住所有格式。5.4 坐标系误判导致点位偏移现象点画在地图上整体偏移或者落进海里、落在错误半球。 原因把投影坐标当经纬度用或者数据本身是经纬度但被误当成投影坐标做了to_crs。 解决用total_bounds看数值量级经纬度在百位以内、投影坐标在几十万到几百万。确认是经纬度后再set_crs(EPSG:4326)补上定义需要算距离再to_crs到 UTM。顺序不能反先set_crs声明原始坐标系再to_crs转换。5.5 轨迹线乱窜排序键不完整现象生成的轨迹线像一团乱麻不同台风的点被连到一起。 原因只按时间排序没按编号分组或者编号字段混入空格、类型不一致导致同一台风被拆组。 解决排序键固定为[typhoon_id, dt]排序前对编号做astype(str).str.strip()清洗。生成轨迹后打印track_gdf的条数与groupby的组数对比数量对不上就说明分组有问题。提示这五条里最容易被忽略的是 5.4因为坐标系错的时候图还是能画出来只是位置不对不核对total_bounds很难发现。6. 进阶验证用时间窗口和缓冲区粗判登陆轨迹连完之后我一般会做两件事验证数据质量一是看每个台风的观测时次间隔是否规整二是用陆地范围粗判登陆点。台风最佳路径数据通常是 6 小时一个时次个别加密到 3 小时间隔突然变大说明该时段缺测直接拿去算频次会低估。# 按台风分组算相邻时次间隔小时 clean clean.sort_values([typhoon_id, dt]) clean[gap_h] clean.groupby(typhoon_id)[dt].diff().dt.total_seconds() / 3600 print(clean[gap_h].describe()) # 间隔明显大于 6 小时的行标记为可疑缺测 suspect clean[clean[gap_h] 12] print(跨度过大的记录数:, len(suspect))diff()在每组内做相邻时间差total_seconds()/3600换算成小时。正常间隔应集中在 3 或 6 小时附近describe()的max如果出现几十上百小时说明某些台风中间断档这些台风做强度演变分析时要单独剔除或插值。这一步不需要额外数据纯靠时间列就能做是性价比最高的完整性检查。登陆粗判用陆地多边形和轨迹点做空间关系判断。手头没有陆地边界数据时可以用一个近似的陆地矩形或已有点集做粗略过滤有精确边界时用within或intersects判断哪些点落在陆地内再按时间取第一个落陆点作为登陆时刻。# 假设已有陆地边界 land_gdf任意来源的陆地多边形 # 没有时先用一个大致的经纬度矩形近似仅用于粗筛 from shapely.geometry import box land_approx box(108, 18, 125, 42) # 粗略陆地范围按需调整 landfall clean[clean.within(land_approx)] first_land (landfall.sort_values(dt) .groupby(typhoon_id) .first()[[dt, geometry]]) print(粗判登陆台风数:, len(first_land))within(land_approx)用近似矩形只能粗筛会把近海但未登陆的点也算进来所以结果只适合做量级核对不适合直接当登陆统计口径。要精确判断换成真实海岸线或陆地多边形并在投影坐标下做intersects避免经纬度直接做拓扑判断时的精度问题。groupby().first()取每个台风最早落陆的那条记录sort_values(dt)保证first取的是时间最早而非行序最早——这一点顺序反了登陆时刻就错了。我踩过最狠的一次是拿一份没做时间排序的轨迹直接groupby().first()结果登陆时刻取到了台风消散后的最后一个点统计出来的登陆月份整体后移。从那以后凡是涉及“首次”“最早”的聚合我都先在分组内显式sort_values再取不依赖原始行序。数据本身不会提醒你顺序错了只有交叉核对才会暴露。希望这些步骤能帮你少走几段弯路。本文还有配套的精品资源点击获取