尧图网页设计 
简介面向基站数据分析与区域经济研究场景提供2006—2024年中国内地、香港、澳门、台湾地区的手机基站数据整理包。该数据源自OpenCelliD覆盖全国约183.6万条样本时间跨度近20年字段包括网络类型、网络代数、移动国家/地区码、区域代码、小区标识、经纬度、覆盖范围、测量样本数、平均信号强度、创建及更新日期等可支撑新型信息基础设施与城市韧性、信号覆盖评估、基站密度分析等研究。资源包为单个docx文档大小仅50KB文档内含百度网盘链接与提取码并附有相关参考文献方便读者直接获取300多MB的完整数据文件。目前已有164人学习该资源适合高年级本科生、研究生以及从事空间统计或数字经济研究的科研人员。通过这份材料读者可大幅缩减原始数据搜集与字段梳理时间快速进入建模分析环节文档所附参考文献还能帮助规范引用提升研究效率。1. 中国内地及港澳台手机基站数据300MB 里到底能挖出什么拿到一份覆盖中国内地及港澳台、时间跨度为 2006-2024 年的手机基站数据包含经纬度、覆盖范围和平均信号强度体积约 300MB你会先做什么我一般不会直接画散点图因为散点图只会把“有一堆基站”这件事可视化出来而基站规划要的覆盖率、LBS 要的地理围栏、优化要的弱覆盖区域都需要把这 300MB 转成可聚合、可叠加、可回溯的结构化资产。这套数据的价值在于同时具备时间和空间两个维度时间轴跨越 2G 到 5G 的整个建设周期空间上则覆盖了内地、香港、澳门、台湾多个坐标体系混杂的区域。做网络优化的人可以看平均信号强度的月度走势做地图产品的人可以把经纬度落到网格上算覆盖盲区做数据治理的人则要先解决格式、编码和坐标偏移三件事。300MB 看起来不大处理起来并不轻松。文件里可能混着 GBK 和 UTF-8 编码经纬度可能是 WGS-84 也可能已经做过 GCJ-02 偏移信号强度字段还会出现 -1、0、-120 这类需要分情况讨论的异常值。下面这套处理流程是我拿到类似数据后最常用的一条路径。2. 解析 300MB 手机基站数据字段识别、编码检查与块式清洗拿到 2006-2024 年的手机基站数据后第一件事不是写 PySpark而是先确认文件到底长什么样。300MB 的体量决定了单机 pandas 完全够用但格式干扰会浪费大量时间。基站数据以 CSV 交付最常见偶尔会碰到制表符分隔、GBK 编码或头部带 BOM 的情况这三种问题都会让read_csv首战失利。2.1 先用三行 shell 命令确认文件类型和字符集在打开文件之前我会先看文件头部字节和编码head -c 2000 cell_tower_2006_2024.csv file -bi cell_tower_2006_2024.csv wc -l cell_tower_2006_2024.csvhead -c 2000显示前 2000 字节几秒内能看到列名和分隔符file -bi输出 MIME 类型和 charset如果显示charsetiso-8859-1不要慌它通常意味着文件里混有 GBK 编码的中文区域名wc -l给出总行数用来估算 chunk 大小。如果拿到的是 gzip 压缩包把head换成zcat cell_tower_2006_2024.csv.gz | head -c 2000其他命令同理。这步的价值在于“先建 schema 再读数据”。对于跨 19 年的基站数据常见字段包括小区标识、经纬度、覆盖半径、平均信号强度和记录时间但不同年份导出的列名可能完全不同。先看头文件再决定下面对字段做 rename 还是直接按位置读。我在处理这类数据时见过的典型字段如下具体名称可能不同但基本可以映射字段示例类型说明cell_id686233int小区编号跨年可能复用lac/tac4201int位置区码用于定位关联longitude113.2685float经度GCJ-02 或 WGS-84 需要确认latitude23.1342float纬度avg_signal_dbm-87.5float平均信号强度单位 dBmcoverage_radius_m1200float覆盖半径估计值record_date2018-06-01date数据日期2.2 用 pandas 按块读取并统一数据类型的标准写法字段确认后我不建议一次性pd.read_csv读全量。300MB 虽然不大但配上 object 类型的中文字段后内存会翻两三倍分块读取是更稳妥的做法import pandas as pd CHUNK 500_000 dtype_dict { cell_id: int64, lac: int32, tac: int32, longitude: float32, latitude: float32, avg_signal_dbm: float32, coverage_radius_m: float32, record_date: str, } parts [] for chunk in pd.read_csv( cell_tower_2006_2024.csv, dtypedtype_dict, chunksizeCHUNK, usecolsdtype_dict.keys(), parse_dates[record_date], na_values[, NULL, \\N], encodingutf-8, ): chunk[year_month] chunk[record_date].dt.to_period(M) parts.append(chunk) df pd.concat(parts, ignore_indexTrue) df.to_parquet(cell_tower_clean.parquet, indexFalse)这里chunksize500_000表示每次读 50 万行usecols只加载需要的列dtype_dict把经纬度和信号强度降为 float32内存占用能减少一半。parse_dates[record_date]在分块读时同样生效但要注意如果原文件日期格式不统一会出现整块解析失败此时可以把record_date先读成 str再用pd.to_datetime(..., errorscoerce)单独处理。na_values参数很关键。基站数据导出时经常用\N表示空值如果不指定这些值会变成字符串而不是 NaN后续聚合会报错。chunk 读完后用pd.concat合并再一次性写成 parquet。Parquet 的好处是保留了列类型后续按年月过滤时不需要重新解析。2.3 清洗三条硬规则坐标边界、信号值域和重复记录清洗阶段我用三条硬规则过一遍lat_ok (df[latitude] 3.0) (df[latitude] 54.0) lon_ok (df[longitude] 73.0) (df[longitude] 136.0) signal_ok (df[avg_signal_dbm] -140) (df[avg_signal_dbm] -30) df df[lat_ok lon_ok signal_ok].drop_duplicates( subset[record_date, cell_id, longitude, latitude] ) print(df.groupby(year_month).size().tail())经度范围 73-136、纬度范围 3-54 覆盖中国内地、香港、澳门和台湾的全部实际区域超出这个范围的记录只可能是漂移点或错误坐标。信号强度-30到-140是一个合理的 dBm 区间-30以上的值通常是施工测试产生的假数据-140以下则可能是无服务时的占位值。需要提醒的是不同网络制式的信号范围并不一致比如 LTE 的 RSRP 很少高于-70但如果数据源统一给出的是“平均信号强度”字段先用这个宽区间过滤之后按制式细分才是正确顺序。drop_duplicates的子集要同时包含日期、小区和坐标。同一个小区在一天内多次上报时基站位置不变但平均信号强度会变化保留第一条只适合做覆盖面积计算如果你想统计信号变化应该先按这一对键聚合取均值或中位数。数据清洗完成后2006-2024 的时间跨度就能转化为year_month轴上的一个连续序列下一章开始处理坐标本身。3. 经纬度坐标转换与覆盖范围网格化让基站数据能直接叠加到地图上手机基站数据里的经纬度不一定是真实 GPS 坐标。在中国内地商用地图坐标系中公开地图和定位接口普遍使用 GCJ-02 加密坐标而基站采集设备输出的往往是 WGS-84。如果不做判断直接把经纬度画到高德底图上会出现几十到几百米的整体偏移基站落在道路对面甚至河对岸。3.1 先判断经纬度是 WGS-84 还是加偏后的 GCJ-02判断方法很朴素取一个你确定位置的基站或地标对比数据里的坐标和地图显示的位置。如果所有点都朝着同一个方向偏移说明整份数据可能是 WGS-84需要转换成 GCJ-02 才能和地图匹配如果点与地图吻合说明数据已经转过了继续转换等于二次加偏。提示如果原始数据里的列名已经带有 gcj02 字样或者你确认数据来自内地商用地图 SDK不要执行 WGS-84 转换。常用的wgs84_to_gcj02转换函数是这样一个实现项目里可以直接拷贝import math def wgs84_to_gcj02(lng, lat): a 6378245.0 ee 0.006693421622965943 if 73 lng or lng 136 or 3 lat or lat 54: return lng, lat d_lat transform_lat(lng - 105.0, lat - 35.0) d_lng transform_lng(lng - 105.0, lat - 35.0) rad_lat lat / 180.0 * math.pi magic math.sin(rad_lat) magic 1 - ee * magic * magic sqrt_magic math.sqrt(magic) d_lat (d_lat * 180.0) / ((a * (1 - ee)) / (magic * sqrt_magic) * math.pi) d_lng (d_lng * 180.0) / (a / sqrt_magic * math.cos(rad_lat) * math.pi) return lng d_lng, lat d_lat def transform_lng(x, y): return 300.0 x 2.0 * y 0.1 * x * x 0.1 * x * y 0.1 * math.sqrt(abs(x)) (20.0 * math.sin(6.0 * x * math.pi) 20.0 * math.sin(2.0 * x * math.pi)) * 2.0 / 3.0 def transform_lat(x, y): return -100.0 2.0 * x 3.0 * y 0.2 * y * y 0.1 * x * y 0.2 * math.sqrt(abs(x)) (20.0 * math.sin(6.0 * x * math.pi) 20.0 * math.sin(2.0 * x * math.pi)) * 2.0 / 3.0 (20.0 * math.sin(y * math.pi) 40.0 * math.sin(y / 3.0 * math.pi)) * 2.0 / 3.0 (160.0 * math.sin(y / 12.0 * math.pi) 320.0 * math.sin(y * math.pi / 30.0)) * 2.0 / 3.0这个公式并不是数学上的精密投影而是对偏移规律的经验拟合。a和ee是坐标转换中固定使用的椭球参数transform_lng和transform_lat负责计算在某个位置的偏移量。实际项目中直接把lng, lat传给函数即可香港、澳门、台湾在内的全境坐标都会落在这个转换的适用范围里。必须注意如果数据里已经有一部分点转过了一部分没转转换就不能全表执行。我会先按地区抽样 500 个点分别用原始值和转换值去叠地图统计偏移方向的方差来判断是否混用。3.2 用 H3 六边形网格把点坐标变成覆盖范围覆盖范围不能直接用圆来画。基站信号受地形、建筑物和天线方位影响实际覆盖范围是不规则多边形但对 300MB 数据做精确建模并不现实。常见做法是用 Uber H3 六边形网格对二维空间做离散化把经纬度转成网格编号之后所有覆盖范围和信号强度统计都基于网格进行。import h3 import pandas as pd df[h3_09] df.apply( lambda r: h3.latlng_to_cell(r[latitude], r[longitude], 9), axis1 ) agg df.groupby([year_month, h3_09]).agg( tower_cnt(cell_id, nunique), avg_signal(avg_signal_dbm, mean), med_radius(coverage_radius_m, median), ).reset_index() agg[wkt] agg[h3_09].map(h3.cell_to_boundary)latlng_to_cell的第三个参数是resolution分辨率越高网格越密。resolution 9 的六边形边长约 400 米级别适合全国宏观覆盖评估要细到一个街道或者大型园区的盲区可以提升到 10 或 11。cell_to_boundary返回多边形的经纬度拐点组成 WKT 后可以直接写入 PostGIS 或 GeoJSON。在 19 年数据上逐行调用apply会比较慢。可以先用pandas.factorize把经纬度组合成整数索引再对唯一坐标对调用 H3 转换之后map回 DataFrame当数据里同一个基站在多年月度重复出现 100 万次时这能省下大部分时间。3.3 批量补齐地址与海拔给基站坐标增加可用属性基站数据只有经纬度时很多分析没法下探比如“弱覆盖区域是在山谷里还是平地上”。这时需要给每个基站批量补地址或高程。大多数公开 API 都有并发限制我会用asyncio做一个可控并发的巡查脚本import asyncio import aiohttp async def fetch_elevation(session, lng, lat): url fhttps://api.open-elevation.com/api/v1/lookup?locations{lat},{lng} async with session.get(url, timeoutaiohttp.ClientTimeout(total10)) as resp: data await resp.json() return data[results][0][elevation] async def batch_elevation(coords, limit20): connector aiohttp.TCPConnector(limitlimit) async with aiohttp.ClientSession(connectorconnector) as session: return await asyncio.gather(*[fetch_elevation(session, lng, lat) for lng, lat in coords])TCPConnector(limit20)控制最多同时建立 20 个连接避免把免费服务打爆。高程数据返回的是 WGS-84 坐标系下的海拔和基站坐标属于同一套基准不需要再做转换。如果数据集中在香港、澳门等山地城市高程对覆盖半径的修正效果非常明显如果只是做全国宏观汇总也不建议跳过因为海拔每上升 100 米基站覆盖半径的变化可能超过 20%。逆地址补全也一样先用高德逆地理编码接口传入已经转好的 GCJ-02 坐标返回省市区和街道名。千万不要把 WGS-84 坐标直接传给内地地图接口否则返回的地址可能与实际位置相差一条街区。4. 平均信号强度聚合与覆盖范围计算从 dBm 统计到热力图输出信号强度是整个数据集里最有业务价值的字段但它也是最容易被误导的字段。先理解口径再动手聚合。4.1 先按网络制式分开统计平均信号强度2006-2024 年横跨 2G、3G、4G、5G不同制式上报的“信号强度”并非同一个物理量。GSM 上报的是 RSSILTE 上报的是 RSRP5G 则是 SS-RSRP三者的正常区间和业务含义都不同。如果数据里有network_type字段算平均信号强度前一定要分组制式指标典型区间覆盖半径参考GSMRSSI-110 ~ -55 dBm几百米到几公里WCDMARSCP-120 ~ -70 dBm城区 300-800 米LTERSRP-140 ~ -75 dBm城区 300-600 米5G NRSS-RSRP-140 ~ -75 dBm城区 150-400 米如果数据源没有制式字段也可以根据时间和经纬度做一个间接判断。2006-2012 年出现的基站大概率是 2G 或 3G 站点2015 年以后新出现在城市高密度区域的站点更可能是 LTE。简单混合所有制式算统一均值会把真正弱覆盖的 LTE 站点和信号良好的 GSM 站点互相抵消。分组统计的代码很直接year_signal df.groupby([year, network_type]).agg( mean_signal(avg_signal_dbm, mean), median_signal(avg_signal_dbm, median), cell_count(cell_id, nunique), ).reset_index() print(year_signal[year_signal[year] 2015].head(20))这里用median_signal而不是只信mean_signal因为基站数据的信号强度分布通常是左偏的极少数远距离弱信号会把平均值拉低。cell_count用于观察 19 年中基站规模的时间拐点如果某一年数量出现断崖式下跌要回去查数据采集覆盖率而不是真的发生了大规模退网。4.2 用 scipy 插值把离散基站信号变成连续覆盖面网格化统计后每个六边形格网里有一个平均信号强度但真正需要输出“覆盖范围图”时格网边界太生硬。更顺滑的做法是使用scipy.interpolate.griddata做空间插值把离散点拟合成一个连续的面。import numpy as np from scipy.interpolate import griddata sample df[ (df[year_month] 2023-06) df[avg_signal_dbm].notna() ].copy() lng_grid, lat_grid np.meshgrid( np.linspace(sample[longitude].min(), sample[longitude].max(), 500), np.linspace(sample[latitude].min(), sample[latitude].max(), 400), ) signal_surface griddata( pointssample[[longitude, latitude]].values, valuessample[avg_signal_dbm].values, xi(lng_grid, lat_grid), methodlinear, fill_value-120, )np.linspace分成 500 x 400 20 万个网格点这个规模在任何笔记本上都能快速出结果。methodlinear速度快但会出现以插值点为中心的三角形折痕想要更平滑就用cubic但一旦点数超过几万cubic 的内存开销会迅速上升。fill_value-120是给无法插值的空白区域一个弱信号兜底值后续做阈值提取时不会误判为盲区外的空洞。全国范围一次插值通常会吃掉好几个 GB 内存所以我一般在上面这段代码之前先按城市或区域过滤掉一部分网格例如只保留香港、澳门城市群或者用 H3 resolution 10 的网格聚合后先降采样再插值。19 年数据逐月插值不现实正确顺序是先按年筛选或者只针对需要出报告的那几个月。4.3 用阈值提取弱覆盖盲区并导出 GeoJSON插值完成后覆盖率计算的本质是一次阈值判断。例如把平均信号强度低于-95 dBm的区域视为弱覆盖低于-105 dBm视为覆盖盲区。提取结果可以用 GeoJSON 直接交给 GIS 部门。import geopandas as gpd from shapely.geometry import Point threshold -95 weak_mask signal_surface threshold weak_points gpd.GeoDataFrame( { longitude: lng_grid[weak_mask], latitude: lat_grid[weak_mask], signal_dbm: signal_surface[weak_mask], }, geometry[ Point(x, y) for x, y in zip(lng_grid[weak_mask], lat_grid[weak_mask]) ], crsEPSG:4326, ) weak_points.to_file(weak_cover_2023_06.geojson, driverGeoJSON)lng_grid[weak_mask]拿到弱覆盖点的经纬度索引再构造成 shapely 的 Point 对象。crsEPSG:4326表示这是 WGS-84 经纬度坐标和前面 H3 网格用的是同一个空间基准。注意如果原始点只集中在市区插值会把郊区大片无数据区域填充为弱覆盖造成“盲区蔓延”的假象。我建议先对样本点做一遍凸包只保留凸包内的插值结果凸包外统一标记为“无采样数据”而不是归为弱覆盖。数据覆盖 2006-2024 年跨度大时间越早基站越稀疏盲区判断越要保守。5. 让 300MB 基站数据可长期复用Parquet 分区、坐标基准留存与栅格缓存前四章跑完后你会得到一份清洗过的 Parquet、一张坐标转换函数表、一张信号强度热力图。这还不够我通常会在项目收尾前做三件小事把 300MB 原始数据的复用成本再降一截。第一件小事是 Parquet 分区。把cell_tower_clean.parquet重写成年份分区df[year] df[record_date].dt.year df.to_parquet(cell_tower_pq/, partition_cols[year], compressionzstd)分区后查询某一年数据时DuckDB 或 Spark 可以跳过其他年份的文件。19 年数据被拆成 19 个目录单次查询可能只读其中两三年。zstd压缩比高300MB 的 CSV 转换后通常只有几十 MB配合 DuckDB 查询可以做到秒级响应import duckdb res duckdb.sql( SELECT year, count(*) AS cnt FROM cell_tower_pq/*.parquet WHERE longitude BETWEEN 113 AND 115 GROUP BY year ORDER BY year ).df()第二件小事是保留双坐标列。清洗阶段不要只存gcj02_lng/gcj02_lat要把原始 WGS-84 经纬度也留在同一张表里。原因是不同下游系统要求不同地图可视化用 GCJ-02基站仿真和传播模型通常用 WGS-84而原始坐标一旦丢弃很难从加偏后的坐标反向恢复。我通常会在 parquet 里同时放wgs84_lng、wgs84_lat、gcj02_lng、gcj02_lat四列宁可多占一点存储也不要让下次对接时重新做抽样判断。第三件小事是栅格缓存。scipy griddata插值在大范围数据上很耗时热力图生成后直接写成 GeoTIFFimport rasterio from rasterio.transform import from_origin with rasterio.open( signal_2023_06.tif, w, driverGTiff, heightsignal_surface.shape[0], widthsignal_surface.shape[1], count1, dtypefloat32, crsEPSG:4326, transformfrom_origin(lng_grid[0][0], lat_grid[0][0], 0.01, 0.01), ) as dst: dst.write(signal_surface.astype(float32), 1)from_origin(west, north, pixel_size_x, pixel_size_y)是左上角坐标加像素尺寸的标准写法0.01 度约等于 1 公里。缓存文件生成后后续覆盖率判断、盲区提取都直接读 TIFF不再重跑插值。把年份分区、双坐标列和 GeoTIFF 缓存放在同一个数据目录里300MB 的原始 CSV 就不再是唯一真相后续每次分析都从 Parquet 和栅格文件出发省掉重复解码和插值的时间。本文还有配套的精品资源点击获取 手机基站数据清洗与覆盖分析实战：从300MB到可复用资产
简介面向基站数据分析与区域经济研究场景提供2006—2024年中国内地、香港、澳门、台湾地区的手机基站数据整理包。该数据源自OpenCelliD覆盖全国约183.6万条样本时间跨度近20年字段包括网络类型、网络代数、移动国家/地区码、区域代码、小区标识、经纬度、覆盖范围、测量样本数、平均信号强度、创建及更新日期等可支撑新型信息基础设施与城市韧性、信号覆盖评估、基站密度分析等研究。资源包为单个docx文档大小仅50KB文档内含百度网盘链接与提取码并附有相关参考文献方便读者直接获取300多MB的完整数据文件。目前已有164人学习该资源适合高年级本科生、研究生以及从事空间统计或数字经济研究的科研人员。通过这份材料读者可大幅缩减原始数据搜集与字段梳理时间快速进入建模分析环节文档所附参考文献还能帮助规范引用提升研究效率。1. 中国内地及港澳台手机基站数据300MB 里到底能挖出什么拿到一份覆盖中国内地及港澳台、时间跨度为 2006-2024 年的手机基站数据包含经纬度、覆盖范围和平均信号强度体积约 300MB你会先做什么我一般不会直接画散点图因为散点图只会把“有一堆基站”这件事可视化出来而基站规划要的覆盖率、LBS 要的地理围栏、优化要的弱覆盖区域都需要把这 300MB 转成可聚合、可叠加、可回溯的结构化资产。这套数据的价值在于同时具备时间和空间两个维度时间轴跨越 2G 到 5G 的整个建设周期空间上则覆盖了内地、香港、澳门、台湾多个坐标体系混杂的区域。做网络优化的人可以看平均信号强度的月度走势做地图产品的人可以把经纬度落到网格上算覆盖盲区做数据治理的人则要先解决格式、编码和坐标偏移三件事。300MB 看起来不大处理起来并不轻松。文件里可能混着 GBK 和 UTF-8 编码经纬度可能是 WGS-84 也可能已经做过 GCJ-02 偏移信号强度字段还会出现 -1、0、-120 这类需要分情况讨论的异常值。下面这套处理流程是我拿到类似数据后最常用的一条路径。2. 解析 300MB 手机基站数据字段识别、编码检查与块式清洗拿到 2006-2024 年的手机基站数据后第一件事不是写 PySpark而是先确认文件到底长什么样。300MB 的体量决定了单机 pandas 完全够用但格式干扰会浪费大量时间。基站数据以 CSV 交付最常见偶尔会碰到制表符分隔、GBK 编码或头部带 BOM 的情况这三种问题都会让read_csv首战失利。2.1 先用三行 shell 命令确认文件类型和字符集在打开文件之前我会先看文件头部字节和编码head -c 2000 cell_tower_2006_2024.csv file -bi cell_tower_2006_2024.csv wc -l cell_tower_2006_2024.csvhead -c 2000显示前 2000 字节几秒内能看到列名和分隔符file -bi输出 MIME 类型和 charset如果显示charsetiso-8859-1不要慌它通常意味着文件里混有 GBK 编码的中文区域名wc -l给出总行数用来估算 chunk 大小。如果拿到的是 gzip 压缩包把head换成zcat cell_tower_2006_2024.csv.gz | head -c 2000其他命令同理。这步的价值在于“先建 schema 再读数据”。对于跨 19 年的基站数据常见字段包括小区标识、经纬度、覆盖半径、平均信号强度和记录时间但不同年份导出的列名可能完全不同。先看头文件再决定下面对字段做 rename 还是直接按位置读。我在处理这类数据时见过的典型字段如下具体名称可能不同但基本可以映射字段示例类型说明cell_id686233int小区编号跨年可能复用lac/tac4201int位置区码用于定位关联longitude113.2685float经度GCJ-02 或 WGS-84 需要确认latitude23.1342float纬度avg_signal_dbm-87.5float平均信号强度单位 dBmcoverage_radius_m1200float覆盖半径估计值record_date2018-06-01date数据日期2.2 用 pandas 按块读取并统一数据类型的标准写法字段确认后我不建议一次性pd.read_csv读全量。300MB 虽然不大但配上 object 类型的中文字段后内存会翻两三倍分块读取是更稳妥的做法import pandas as pd CHUNK 500_000 dtype_dict { cell_id: int64, lac: int32, tac: int32, longitude: float32, latitude: float32, avg_signal_dbm: float32, coverage_radius_m: float32, record_date: str, } parts [] for chunk in pd.read_csv( cell_tower_2006_2024.csv, dtypedtype_dict, chunksizeCHUNK, usecolsdtype_dict.keys(), parse_dates[record_date], na_values[, NULL, \\N], encodingutf-8, ): chunk[year_month] chunk[record_date].dt.to_period(M) parts.append(chunk) df pd.concat(parts, ignore_indexTrue) df.to_parquet(cell_tower_clean.parquet, indexFalse)这里chunksize500_000表示每次读 50 万行usecols只加载需要的列dtype_dict把经纬度和信号强度降为 float32内存占用能减少一半。parse_dates[record_date]在分块读时同样生效但要注意如果原文件日期格式不统一会出现整块解析失败此时可以把record_date先读成 str再用pd.to_datetime(..., errorscoerce)单独处理。na_values参数很关键。基站数据导出时经常用\N表示空值如果不指定这些值会变成字符串而不是 NaN后续聚合会报错。chunk 读完后用pd.concat合并再一次性写成 parquet。Parquet 的好处是保留了列类型后续按年月过滤时不需要重新解析。2.3 清洗三条硬规则坐标边界、信号值域和重复记录清洗阶段我用三条硬规则过一遍lat_ok (df[latitude] 3.0) (df[latitude] 54.0) lon_ok (df[longitude] 73.0) (df[longitude] 136.0) signal_ok (df[avg_signal_dbm] -140) (df[avg_signal_dbm] -30) df df[lat_ok lon_ok signal_ok].drop_duplicates( subset[record_date, cell_id, longitude, latitude] ) print(df.groupby(year_month).size().tail())经度范围 73-136、纬度范围 3-54 覆盖中国内地、香港、澳门和台湾的全部实际区域超出这个范围的记录只可能是漂移点或错误坐标。信号强度-30到-140是一个合理的 dBm 区间-30以上的值通常是施工测试产生的假数据-140以下则可能是无服务时的占位值。需要提醒的是不同网络制式的信号范围并不一致比如 LTE 的 RSRP 很少高于-70但如果数据源统一给出的是“平均信号强度”字段先用这个宽区间过滤之后按制式细分才是正确顺序。drop_duplicates的子集要同时包含日期、小区和坐标。同一个小区在一天内多次上报时基站位置不变但平均信号强度会变化保留第一条只适合做覆盖面积计算如果你想统计信号变化应该先按这一对键聚合取均值或中位数。数据清洗完成后2006-2024 的时间跨度就能转化为year_month轴上的一个连续序列下一章开始处理坐标本身。3. 经纬度坐标转换与覆盖范围网格化让基站数据能直接叠加到地图上手机基站数据里的经纬度不一定是真实 GPS 坐标。在中国内地商用地图坐标系中公开地图和定位接口普遍使用 GCJ-02 加密坐标而基站采集设备输出的往往是 WGS-84。如果不做判断直接把经纬度画到高德底图上会出现几十到几百米的整体偏移基站落在道路对面甚至河对岸。3.1 先判断经纬度是 WGS-84 还是加偏后的 GCJ-02判断方法很朴素取一个你确定位置的基站或地标对比数据里的坐标和地图显示的位置。如果所有点都朝着同一个方向偏移说明整份数据可能是 WGS-84需要转换成 GCJ-02 才能和地图匹配如果点与地图吻合说明数据已经转过了继续转换等于二次加偏。提示如果原始数据里的列名已经带有 gcj02 字样或者你确认数据来自内地商用地图 SDK不要执行 WGS-84 转换。常用的wgs84_to_gcj02转换函数是这样一个实现项目里可以直接拷贝import math def wgs84_to_gcj02(lng, lat): a 6378245.0 ee 0.006693421622965943 if 73 lng or lng 136 or 3 lat or lat 54: return lng, lat d_lat transform_lat(lng - 105.0, lat - 35.0) d_lng transform_lng(lng - 105.0, lat - 35.0) rad_lat lat / 180.0 * math.pi magic math.sin(rad_lat) magic 1 - ee * magic * magic sqrt_magic math.sqrt(magic) d_lat (d_lat * 180.0) / ((a * (1 - ee)) / (magic * sqrt_magic) * math.pi) d_lng (d_lng * 180.0) / (a / sqrt_magic * math.cos(rad_lat) * math.pi) return lng d_lng, lat d_lat def transform_lng(x, y): return 300.0 x 2.0 * y 0.1 * x * x 0.1 * x * y 0.1 * math.sqrt(abs(x)) (20.0 * math.sin(6.0 * x * math.pi) 20.0 * math.sin(2.0 * x * math.pi)) * 2.0 / 3.0 def transform_lat(x, y): return -100.0 2.0 * x 3.0 * y 0.2 * y * y 0.1 * x * y 0.2 * math.sqrt(abs(x)) (20.0 * math.sin(6.0 * x * math.pi) 20.0 * math.sin(2.0 * x * math.pi)) * 2.0 / 3.0 (20.0 * math.sin(y * math.pi) 40.0 * math.sin(y / 3.0 * math.pi)) * 2.0 / 3.0 (160.0 * math.sin(y / 12.0 * math.pi) 320.0 * math.sin(y * math.pi / 30.0)) * 2.0 / 3.0这个公式并不是数学上的精密投影而是对偏移规律的经验拟合。a和ee是坐标转换中固定使用的椭球参数transform_lng和transform_lat负责计算在某个位置的偏移量。实际项目中直接把lng, lat传给函数即可香港、澳门、台湾在内的全境坐标都会落在这个转换的适用范围里。必须注意如果数据里已经有一部分点转过了一部分没转转换就不能全表执行。我会先按地区抽样 500 个点分别用原始值和转换值去叠地图统计偏移方向的方差来判断是否混用。3.2 用 H3 六边形网格把点坐标变成覆盖范围覆盖范围不能直接用圆来画。基站信号受地形、建筑物和天线方位影响实际覆盖范围是不规则多边形但对 300MB 数据做精确建模并不现实。常见做法是用 Uber H3 六边形网格对二维空间做离散化把经纬度转成网格编号之后所有覆盖范围和信号强度统计都基于网格进行。import h3 import pandas as pd df[h3_09] df.apply( lambda r: h3.latlng_to_cell(r[latitude], r[longitude], 9), axis1 ) agg df.groupby([year_month, h3_09]).agg( tower_cnt(cell_id, nunique), avg_signal(avg_signal_dbm, mean), med_radius(coverage_radius_m, median), ).reset_index() agg[wkt] agg[h3_09].map(h3.cell_to_boundary)latlng_to_cell的第三个参数是resolution分辨率越高网格越密。resolution 9 的六边形边长约 400 米级别适合全国宏观覆盖评估要细到一个街道或者大型园区的盲区可以提升到 10 或 11。cell_to_boundary返回多边形的经纬度拐点组成 WKT 后可以直接写入 PostGIS 或 GeoJSON。在 19 年数据上逐行调用apply会比较慢。可以先用pandas.factorize把经纬度组合成整数索引再对唯一坐标对调用 H3 转换之后map回 DataFrame当数据里同一个基站在多年月度重复出现 100 万次时这能省下大部分时间。3.3 批量补齐地址与海拔给基站坐标增加可用属性基站数据只有经纬度时很多分析没法下探比如“弱覆盖区域是在山谷里还是平地上”。这时需要给每个基站批量补地址或高程。大多数公开 API 都有并发限制我会用asyncio做一个可控并发的巡查脚本import asyncio import aiohttp async def fetch_elevation(session, lng, lat): url fhttps://api.open-elevation.com/api/v1/lookup?locations{lat},{lng} async with session.get(url, timeoutaiohttp.ClientTimeout(total10)) as resp: data await resp.json() return data[results][0][elevation] async def batch_elevation(coords, limit20): connector aiohttp.TCPConnector(limitlimit) async with aiohttp.ClientSession(connectorconnector) as session: return await asyncio.gather(*[fetch_elevation(session, lng, lat) for lng, lat in coords])TCPConnector(limit20)控制最多同时建立 20 个连接避免把免费服务打爆。高程数据返回的是 WGS-84 坐标系下的海拔和基站坐标属于同一套基准不需要再做转换。如果数据集中在香港、澳门等山地城市高程对覆盖半径的修正效果非常明显如果只是做全国宏观汇总也不建议跳过因为海拔每上升 100 米基站覆盖半径的变化可能超过 20%。逆地址补全也一样先用高德逆地理编码接口传入已经转好的 GCJ-02 坐标返回省市区和街道名。千万不要把 WGS-84 坐标直接传给内地地图接口否则返回的地址可能与实际位置相差一条街区。4. 平均信号强度聚合与覆盖范围计算从 dBm 统计到热力图输出信号强度是整个数据集里最有业务价值的字段但它也是最容易被误导的字段。先理解口径再动手聚合。4.1 先按网络制式分开统计平均信号强度2006-2024 年横跨 2G、3G、4G、5G不同制式上报的“信号强度”并非同一个物理量。GSM 上报的是 RSSILTE 上报的是 RSRP5G 则是 SS-RSRP三者的正常区间和业务含义都不同。如果数据里有network_type字段算平均信号强度前一定要分组制式指标典型区间覆盖半径参考GSMRSSI-110 ~ -55 dBm几百米到几公里WCDMARSCP-120 ~ -70 dBm城区 300-800 米LTERSRP-140 ~ -75 dBm城区 300-600 米5G NRSS-RSRP-140 ~ -75 dBm城区 150-400 米如果数据源没有制式字段也可以根据时间和经纬度做一个间接判断。2006-2012 年出现的基站大概率是 2G 或 3G 站点2015 年以后新出现在城市高密度区域的站点更可能是 LTE。简单混合所有制式算统一均值会把真正弱覆盖的 LTE 站点和信号良好的 GSM 站点互相抵消。分组统计的代码很直接year_signal df.groupby([year, network_type]).agg( mean_signal(avg_signal_dbm, mean), median_signal(avg_signal_dbm, median), cell_count(cell_id, nunique), ).reset_index() print(year_signal[year_signal[year] 2015].head(20))这里用median_signal而不是只信mean_signal因为基站数据的信号强度分布通常是左偏的极少数远距离弱信号会把平均值拉低。cell_count用于观察 19 年中基站规模的时间拐点如果某一年数量出现断崖式下跌要回去查数据采集覆盖率而不是真的发生了大规模退网。4.2 用 scipy 插值把离散基站信号变成连续覆盖面网格化统计后每个六边形格网里有一个平均信号强度但真正需要输出“覆盖范围图”时格网边界太生硬。更顺滑的做法是使用scipy.interpolate.griddata做空间插值把离散点拟合成一个连续的面。import numpy as np from scipy.interpolate import griddata sample df[ (df[year_month] 2023-06) df[avg_signal_dbm].notna() ].copy() lng_grid, lat_grid np.meshgrid( np.linspace(sample[longitude].min(), sample[longitude].max(), 500), np.linspace(sample[latitude].min(), sample[latitude].max(), 400), ) signal_surface griddata( pointssample[[longitude, latitude]].values, valuessample[avg_signal_dbm].values, xi(lng_grid, lat_grid), methodlinear, fill_value-120, )np.linspace分成 500 x 400 20 万个网格点这个规模在任何笔记本上都能快速出结果。methodlinear速度快但会出现以插值点为中心的三角形折痕想要更平滑就用cubic但一旦点数超过几万cubic 的内存开销会迅速上升。fill_value-120是给无法插值的空白区域一个弱信号兜底值后续做阈值提取时不会误判为盲区外的空洞。全国范围一次插值通常会吃掉好几个 GB 内存所以我一般在上面这段代码之前先按城市或区域过滤掉一部分网格例如只保留香港、澳门城市群或者用 H3 resolution 10 的网格聚合后先降采样再插值。19 年数据逐月插值不现实正确顺序是先按年筛选或者只针对需要出报告的那几个月。4.3 用阈值提取弱覆盖盲区并导出 GeoJSON插值完成后覆盖率计算的本质是一次阈值判断。例如把平均信号强度低于-95 dBm的区域视为弱覆盖低于-105 dBm视为覆盖盲区。提取结果可以用 GeoJSON 直接交给 GIS 部门。import geopandas as gpd from shapely.geometry import Point threshold -95 weak_mask signal_surface threshold weak_points gpd.GeoDataFrame( { longitude: lng_grid[weak_mask], latitude: lat_grid[weak_mask], signal_dbm: signal_surface[weak_mask], }, geometry[ Point(x, y) for x, y in zip(lng_grid[weak_mask], lat_grid[weak_mask]) ], crsEPSG:4326, ) weak_points.to_file(weak_cover_2023_06.geojson, driverGeoJSON)lng_grid[weak_mask]拿到弱覆盖点的经纬度索引再构造成 shapely 的 Point 对象。crsEPSG:4326表示这是 WGS-84 经纬度坐标和前面 H3 网格用的是同一个空间基准。注意如果原始点只集中在市区插值会把郊区大片无数据区域填充为弱覆盖造成“盲区蔓延”的假象。我建议先对样本点做一遍凸包只保留凸包内的插值结果凸包外统一标记为“无采样数据”而不是归为弱覆盖。数据覆盖 2006-2024 年跨度大时间越早基站越稀疏盲区判断越要保守。5. 让 300MB 基站数据可长期复用Parquet 分区、坐标基准留存与栅格缓存前四章跑完后你会得到一份清洗过的 Parquet、一张坐标转换函数表、一张信号强度热力图。这还不够我通常会在项目收尾前做三件小事把 300MB 原始数据的复用成本再降一截。第一件小事是 Parquet 分区。把cell_tower_clean.parquet重写成年份分区df[year] df[record_date].dt.year df.to_parquet(cell_tower_pq/, partition_cols[year], compressionzstd)分区后查询某一年数据时DuckDB 或 Spark 可以跳过其他年份的文件。19 年数据被拆成 19 个目录单次查询可能只读其中两三年。zstd压缩比高300MB 的 CSV 转换后通常只有几十 MB配合 DuckDB 查询可以做到秒级响应import duckdb res duckdb.sql( SELECT year, count(*) AS cnt FROM cell_tower_pq/*.parquet WHERE longitude BETWEEN 113 AND 115 GROUP BY year ORDER BY year ).df()第二件小事是保留双坐标列。清洗阶段不要只存gcj02_lng/gcj02_lat要把原始 WGS-84 经纬度也留在同一张表里。原因是不同下游系统要求不同地图可视化用 GCJ-02基站仿真和传播模型通常用 WGS-84而原始坐标一旦丢弃很难从加偏后的坐标反向恢复。我通常会在 parquet 里同时放wgs84_lng、wgs84_lat、gcj02_lng、gcj02_lat四列宁可多占一点存储也不要让下次对接时重新做抽样判断。第三件小事是栅格缓存。scipy griddata插值在大范围数据上很耗时热力图生成后直接写成 GeoTIFFimport rasterio from rasterio.transform import from_origin with rasterio.open( signal_2023_06.tif, w, driverGTiff, heightsignal_surface.shape[0], widthsignal_surface.shape[1], count1, dtypefloat32, crsEPSG:4326, transformfrom_origin(lng_grid[0][0], lat_grid[0][0], 0.01, 0.01), ) as dst: dst.write(signal_surface.astype(float32), 1)from_origin(west, north, pixel_size_x, pixel_size_y)是左上角坐标加像素尺寸的标准写法0.01 度约等于 1 公里。缓存文件生成后后续覆盖率判断、盲区提取都直接读 TIFF不再重跑插值。把年份分区、双坐标列和 GeoTIFF 缓存放在同一个数据目录里300MB 的原始 CSV 就不再是唯一真相后续每次分析都从 Parquet 和栅格文件出发省掉重复解码和插值的时间。本文还有配套的精品资源点击获取