拓冰网站定制 
1. 为什么我坚持用 Jupyter Notebook 做数据预处理先交代一下背景。我日常工作中有一大半时间花在“把乱七八糟的数据弄成能用的数据”——这一步有个专业名字叫数据预处理。以前我习惯用 PyCharm 写一个完整的 .py 文件跑完看结果错了一步就得从头看日志改一处变量又要重新跑一遍全流程。直到我尝试把整个流程搬到 Jupyter Notebook 上才发现自己以前浪费了多少时间。如果在看这篇文章的你也是个经常和数据打交道的人我建议你认真看完因为数据预处理场景下交互式执行的优势是常态化的刚需而不是锦上添花。Jupyter Notebook 的核心价值可以概括为三个词分段执行、状态保留、可视化反馈。分段执行意味着我可以把“读数据—清洗—转换—统计—导出”拆成一个一个的单元格跑到哪里断了就在哪里修不用每次从头跑。状态保留则解决了“内存里的东西不丢”的问题你上一步读入的 DataFrame、计算好的掩膜、调整过的参数全部存在内核里后续单元格能直接引用。可视化反馈更是杀手锏一个matplotlib出图放在那里数据长什么样一目了然比对着 10 万行的表格看头部几行要直观得多。这套工作流不仅适合表格数据对栅格数据、遥感影像同样适用。我接下来会用一个 NPP 夜间灯光数据的例子完整走一遍数据预处理的各个环节把安装配置、清洗逻辑、裁剪投影、统计导出、踩坑排查都串起来。你可以把它当成一份直接抄的作业也可以在此基础上替换成你自己的数据。1.1 为什么交互式执行对预处理这么重要说白了数据预处理是个“试错”的过程。你拿到一份数据不知道它的缺失值长什么样、单位是什么、坐标系有没有偏离只能一个个条件去试。脚本语言里一次完整的python script.py执行通常要等几十秒甚至几分钟才能看到结果中间一个字段名写错又要重来。Notebook 则把整个过程变得像“聊天”一样写一个单元格、跑一个单元格立刻能看到输出。我现在的习惯是第一轮先建一个“侦察型 Notebook”专门用来摸数据的底比如 df.info()、df.describe()、df.head()然后是栅格数据的 crs、分辨率、nodata、波段数。这些探索性的判断通常会产生很多临时变量如果在脚本里最后要么删掉要么注释掉非常繁琐在 Notebook 里就无所谓留着也不碍事。第二轮再根据侦察结果写正式的清洗流程。这种“先侦察、后动刀”的思路能减少大量无效代码。还有一个实验室场景的好处就是可以保留实验的“中间状态”。我在处理夜间灯光数据的时候重投影、裁剪这类栅格操作非常耗时经常跑到一半我要去确认某个参数对不对。在 Notebook 里我可以随时暂停改掉参数然后只重跑当前单元格和后续部分前面的数据读入、预处理都不用重新执行。1.2 环境搭建的两种姿势与选择逻辑聊环境总有人问到底是装 Anaconda 还是 Miniconda还是直接用 pip 一个个装。我的看法是如果你不是已经有 Python 基础的老手直接上 Anaconda 最省心。因为 Anaconda 自带 Python 解释器、Jupyter Notebook 和一大堆常用库装完就能用而 Miniconda 只是最小化环境管理器还得手动创建环境和安装包。不要被很多教程里“Anaconda 太笨重”的说法吓到这年头磁盘空间没这么金贵把时间花在安装配置上才是真的不划算。如果一定要用命令行我的推荐是这么一条# 创建独立环境尽量指定 Python 版本避免日后依赖冲突 conda create -n geodata python3.9 -y conda activate geodata # 安装核心数据分析库 conda install -y jupyter pandas numpy matplotlib # 处理栅格、矢量数据的常用库 conda install -y rasterio geopandas rioxarray这里有个细节我明确要求python3.9而不是直接用默认的 Python 3.12 之类。原因很简单rasterio、GDAL 这类底层库的编译依赖非常重虽然官方维护者一直在跟新版本但偶尔还是会遇到某个轮子没跟上最新版 Python 的情况。用 3.9 或者 3.10在 conda-forge 上几乎能找到所有你需要的预编译包少很多麻烦。如果你是 pip 用户也有一条对应的路径pip install notebook pandas numpy matplotlib pip install rasterio geopandas rioxarray不过我个人的经验是GDAL 相关的生态库用 conda 安装的稳定性明显更好。因为 conda 不仅装 Python 包还会统一管理 GDAL、PROJ、GEOS 这些底层 C 库的版本。rasterio 这类库对底层库版本极度敏感pip 安装经常会遇到一个包依赖新版 GDAL、另一个包锁定旧版 GDAL 的冲突。conda 在这方面像个大管家会把所有底层依赖一起协调好。1.3 先装顺手的三件套目录、默认路径和快捷键很多人用 Jupyter Notebook都是打开浏览器新建一个 ipynb然后开始写。但用了一段时间你就会发现不解决三个问题Notebook 用起来会很难受一是文件多了之后目录乱成一团二是不知道文件默认保存到哪里来回切换烦人三是重复操作全靠鼠标效率很低。第一件事是给 Notebook 添加目录。老版本的 Jupyter Notebook比如 6.x可以通过插件实现命令是conda install -c conda-forge jupyter_contrib_nbextensions -y jupyter nbextension enable --py jupyter_contrib_nbextensions --sys-prefix装完重启 Jupyter打开任意一个 ipynb上方菜单栏如果有 Nbextensions 标签进去勾选 Table of Contents (2)页面上就会出现侧边目录。但注意如果你用的是 JupyterLab现在很多新版本的默认界面就是 Lab那么其实不用装任何插件它自带 Table of Contents 面板点左侧那个带层叠图标的面板就能看到当前 notebook 的目录结构。所以如果你在旧版 Notebook 里折腾了半天插件不生效我建议直接切换到 JupyterLab 试试。第二件事修改默认保存路径。默认情况下Jupyter Notebook 会把文件保存在启动时所在目录如果你每次都是从不同的地方启动文件就会散落各处。我自己的做法是把所有 ipynb 统一放在D:/jupyter_workspace下面操作方式分两步# 先生成配置文件 jupyter notebook --generate-config然后用文本编辑器打开生成的文件一般在用户目录下的.jupyter/jupyter_notebook_config.py找到这一行# c.NotebookApp.notebook_dir 改成c.NotebookApp.notebook_dir D:/jupyter_workspace保存之后再启动 Jupyter Notebook你会发现默认工作目录已经变了。这个配置也适用于 JupyterLab配置文件里改成c.ServerApp.notebook_dir也可以看你使用的版本。还有一个小技巧如果你只是想临时指定某个目录不需要改配置可以用jupyter notebook --notebook-dirD:/jupyter_workspace这种方式。第三件事是几个高频快捷键。在 Notebook 里按一次Esc进入命令模式然后B在当前单元格下方新建一个单元格A在当前单元格上方新建一个单元格D D连续按两次 D 删除当前单元格M把当前单元格切换为 MarkdownY把当前单元格切换为代码Shift Enter运行当前单元格并跳到下一个这组快捷键里我最常用的是Shift Enter和M。写数据预处理的文档时我会频繁插入 Markdown 分隔说明用快捷键切换模式比鼠标点击快得多。我的个人心得是数据预处理 Notebook 里Markdown 单元格的注释几乎和代码一样重要你一周后再回头看还能回忆起每一步在处理什么。2. 数据预处理不是写代码是先画流程图很多人处理数据的习惯是拿到数据直接开始写 pandas 操作写一步看一步最后数据处理完了整个流程却说不清楚。这种做法在小数据集上还行一旦数据量上来、步骤变多就很容易出现“处理完发现某个条件判断错了、前面所有清洗都白做”的悲剧。我现在的习惯是动手前先在 Notebook 顶部用 Markdown 画一个流程草图哪怕是几个列表也够用。这一步表面上看耽误了五分钟实际上能帮你节省至少半小时的返工时间。数据预处理说到底是一条流水线原始数据进来经过清洗、变换、合并、采样、校验最终变成可以进入建模或分析的规范数据集。流水线的每一道工序都要想清楚它要解决什么问题。2.1 一次完整预处理的六个步骤如果把数据预处理分成标准步骤大致是这六段数据加载与侦察读入数据了解字段类型、缺失值、分布、坐标系、分辨率等基本属性。数据清洗处理缺失值、重复值、异常值、非法值让数据在“物理层面”变得规范。数据变换重命名、类型转换、归一化、维度调整、重投影、重采样等让数据在“结构层面”变得可用。数据整合多张表拼接、栅格与矢量叠加、时间序列对齐等。数据统计与样本梳理生成描述性统计抽样核对数据合理性。数据导出与文档输出规范文件记录处理参数、版本号和步骤索引。这个框架不只是我的习惯几乎所有成熟的数据团队都有类似的分层。比如说你导出的最终数据出了问题你要能够回溯到它是从第几步开始坏的。如果中间产物都留着排查就很简单如果全都混在一个变量里覆盖来覆盖去就只能从头再来。2.2 把“中间产物规范”定下来我强烈建议在开始预处理之前先定义一个“产物规范”。这里指的是一套固定的命名、格式和组织规则比如原始数据单独放在raw/目录永远不修改修改后的数据放在processed/。中间产物文件名带上步骤号比如01_读入检查.pkl、02_清洗后.csv、03_裁剪.tif、04_重采样.tif。统一用YYYYMMDD作为日期后缀方便按版本回滚。统计结果和清洗后的数据分开存放不要混在一个文件里。这套规范看起来很简单实际操作中却能救你很多次。举个例子我在处理夜间灯光数据时因为中间要经过裁剪、重投影两次栅格运算每步生成的 tif 都很大运行时间也不短。如果我把每一步结果保存成一个中间 tif就算后面的操作出了问题我也只需要从对应的中间结果重新开始而不是把前面的运行又跑一遍。数据预处理有个铁律凡是跑得慢或不好恢复的步骤都值得做中间落盘。2.3 为什么选 NPP 夜间灯光数据做演示案例我做技术选型的时候喜欢挑一个“能覆盖绝大多数问题”的案例数据。NPP 夜间灯光数据恰好就是这么个典型。它是一种栅格影像数据常被用于城市扩展分析、区域社会发展研究等场景优点是数据相对容易获取、覆盖范围广、结构清晰缺点也很明显坐标系需要严格处理、原始数据有大量边界异常值、空间分辨率多种多样、文件体积不小。它的技术难点集中在四个方向一是栅格影像的读写不是寻常的 txt、csv需要 rasterio 这类专业库二是坐标系投影需要特殊处理用错投影直接导致空间偏移三是负值、极值、NoData 这类“脏数据”在影像里比表格里更隐蔽必须通过可视化才能发现四是栅格的统计方式与表格不同需要掩膜、分区统计等专业操作。任何一个玩过遥感或地理数据分析的人都会被这四个问题打磨过一遍。如果你只想快速体验流程可以用任何一张 GeoTIFF 栅格数据替代 NPP 夜间灯光数据步骤完全通用。但如果你想把整个案例完整复现那就跟着我的流程走你会发现这远比单纯处理一个 csv 学到的东西多得多。3. 完整实操NPP 夜间灯光数据预处理从 tif 到统计表下面进入正题。我用一台普通 Windows 笔记本环境是 conda 创建的geodata核心库是rasterio、geopandas、rioxarray。整个案例的目标是把原始 NPP 夜间灯光 GeoTIFF 文件读取、清洗、裁剪、重投影最终输出一个研究区内的像元统计表并生成一张可读的成果图。这里要做一个说明NPP-VIIRS 夜间灯光数据是 NASA/NOAA 合成的月度平均夜间灯光强度产品单位一般为nW/cm2/sr文件格式是 GeoTIFF数据类型有浮点型和整型两种版本。它每个像元的值代表该位置的平均夜间灯光辐亮度值越大通常说明该区域夜间灯光越强。它不像普通 CSV 那样可以直接在 Excel 里打开检查必须借助专业工具。3.1 侦察数据读进来先把“底细”摸清我拿到一份栅格数据从来不会急着写统计代码而是先做一个“读入检查”。这段代码信息量非常大import rasterio import numpy as np import matplotlib.pyplot as plt from rasterio.plot import show # 数据路径自行替换成你自己的文件 tif_path data/raw/VNL_v2_npp_202203_avg_vcmslcfg.tif with rasterio.open(tif_path) as src: print(CRS:, src.crs) print(尺寸(宽,高):, src.width, src.height) print(分辨率:, src.res) print(波段数:, src.count) print(NoData:, src.nodata) print(数据范围:, src.bounds) # 读取第一个波段 band1 src.read(1) print(数据类型:, band1.dtype) print(像元值范围:, np.nanmin(band1), np.nanmax(band1)) # 快速看一眼直方图 plt.figure(figsize(10, 4)) plt.hist(band1[band1 0].ravel(), bins50, range(0, 200)) plt.title(夜间灯光像元值直方图(仅正值)) plt.show()这段代码的输出会告诉我几件事坐标系是什么很多时候 NPP 数据用的是 WGS 84 经纬度坐标分辨率是多少常见的是 15 弧秒大概相当于赤道附近 500 米数据类型的浮点还是整型NoData 用的什么标记。这些信息直接决定了后面的处理参数。很多人在这一步就开始犯错了。有些人直接src.read(1)之后拿着数据就统计结果把 NoData 当成了真实 0 值或者把坐标投影忽略掉直接把经纬度当成平面坐标来算面积和距离。我建议把这一步的打印结果保存下来作为后续处理文档的“数据档案”。3.2 清洗处理负值、NoData 与极值NPP 夜间灯光数据的“脏”主要体现在负值和极端值上。月度合成的 VNL 数据虽然经过了初步处理但依然可能存在少量负值像元这些负值并不是真实灯光而多半是定标、几何校正过程中的噪声产物。我的清洗原则是负值一律先置为零或 NoData不要直接当有效值参与统计。import numpy as np # 重新读取 with rasterio.open(tif_path) as src: band src.read(1) profile src.profile # 查看负值像元数 neg_mask band 0 print(负值像元个数:, np.sum(neg_mask)) print(负值占比: {:.6f}%.format(np.sum(neg_mask) / band.size * 100)) # 清洗负值置为 NoData 或 0 # 这里选择置为0相对于负值更有利于后续统计 band[neg_mask] 0 # 防止极值干扰做一个上限截断 extreme_mask band 200 print(超过200的像元个数:, np.sum(extreme_mask)) # 如果确认是异常值可以截断否则保留 # band[extreme_mask] 200这个上限截断需要谨慎使用。夜间灯光数据在某些高亮度区域比如大型机场、商业中心确实会有很高的值直接截断反而会损失信息。我这里的extreme_mask只是先看一眼数量并不代表一定要处理。实际操作中我会结合研究区域的实际情况判断如果是做全局统计分析且异常像元占比极小我会选择保留如果异常像元挤占直方图导致可视化严重失真我会考虑对数变换或截断。处理异常值没有统一答案核心原则是每一次清洗操作都要记录下来操作原因和参数。NoData 的处理相对标准一些。如果用rioxarray可以这样import rioxarray raster rioxarray.open_rasterio(tif_path, maskedTrue) print(NoData 处理后的数据维度:, raster.shape) # maskedTrue 会把 NoData 位置变成 NaN print(NaN 像元数:, int(raster.isnull().sum()))采用maskedTrue的好处是后续的统计函数可以用nansum、nanmean这类函数直接跳过无效像元不用自己手写掩膜。这一点在处理有大量海域或云覆盖区域的影像时特别重要。3.3 按研究区裁剪并统一坐标系接下来是数据预处理的大头裁剪。也就是只保留你关心区域的像元。如果研究区是一个行政区或一个自定义矩形边界裁剪有两种常见方式按矢量边界裁剪和按经纬度范围裁剪。我以按矢量边界裁剪为例这也是遥感领域最常用的方式。import geopandas as gpd import rasterio from rasterio.mask import mask as rio_mask # 读取研究区矢量边界这里用一份简化面数据 study_area gpd.read_file(data/vector/study_area.shp) # 确保矢量与栅格在同一坐标系 # 如果矢量是 WGS84 经纬度而栅格也是 WGS84则跳过 reproject # 如果不是请用 study_area study_area.to_crs(栅格crs) crs_raster rasterio.open(tif_path).crs study_area study_area.to_crs(crs_raster) # 按边界裁剪 with rasterio.open(tif_path) as src: out_image, out_transform rio_mask(src, study_area.geometry, cropTrue, nodata0) out_meta src.meta.copy() print(裁剪后尺寸:, out_image.shape) # 更新元数据 out_meta.update({ driver: GTiff, height: out_image.shape[1], width: out_image.shape[2], transform: out_transform, nodata: 0 }) # 保存中间结果 with rasterio.open(data/processed/01_裁剪_研究区.tif, w, **out_meta) as dest: dest.write(out_image)这段逻辑里有几个容易踩的坑。第一裁剪之前必须检查矢量与栅格的坐标系是否一致不一致要to_crs转换。第二rio_mask的cropTrue表示对边界进行最小矩形裁剪这样输出影像会比原始影像小很多既能提升后续处理速度也能避免不必要的内存占用。第三nodata0这个参数要和前面清洗策略保持一致如果前面把负值置为 0这里裁掉边界外的部分也用 0 填充整条流水线才一致。如果你只需要一个矩形范围也可以用经纬度范围裁剪代码要简单很多# 用经纬度范围裁剪 min_lon, min_lat 110.0, 30.0 max_lon, max_lat 112.0, 32.0 if raster.rio.crs is None or raster.rio.crs.to_epsg() ! 4326: raster raster.rio.reproject(EPSG:4326) clipped raster.rio.clip_box(minxmin_lon, minymin_lat, maxxmax_lon, maxymax_lat)我这里用的rioxarray它封装了 rasterio很多操作一行就能搞定代码可读性也更好。你如果习惯纯 rasterio那就按第一个示例来两个思路是一样的。3.4 重采样与统计从栅格到表格先同步一个概念重采样就是把不同分辨率的栅格统一到一个目标分辨率。为什么要做这步因为很多时候你要拿夜间灯光数据和别的数据做叠加分析比如人口格网、土地利用影像它们的分辨率可能是 100 米、250 米、1000 米跟 NPP 的 500 米左右不一样。不统一分辨率像元之间无法对齐统计就无从谈起。import rioxarray # 重采样到目标分辨率例如500米 target_resolution 500 clipped rioxarray.open_rasterio(data/processed/01_裁剪_研究区.tif, maskedTrue) # 计算目标网格尺寸 # 先用一个简单的估算方法得到新分辨率下的行列数 width int((clipped.rio.bounds()[2] - clipped.rio.bounds()[0]) / target_resolution) height int((clipped.rio.bounds()[3] - clipped.rio.bounds()[1]) / target_resolution) # 重采样使用平均值方法对夜间灯光这种连续型变量更合适 resampled clipped.rio.reproject( clipped.rio.crs, resolutiontarget_resolution, resamplingrioxarray.enums.Resampling.average ) print(重采样后 shape:, resampled.shape) # 保存重采样结果 resampled.rio.to_raster(data/processed/02_重采样_500m.tif)关于采样算法我简单说一下选择逻辑对连续型数值变量比如夜间灯光辐亮度、气温、降水优先用average平均值重采样因为它不会丢失总辐射信息对分类变量比如土地利用类型优先用nearest最邻近或mode众数因为分类值不能被平均出“疑似新类别”如果是要精确配准、不希望改变像元值的范围可以用bilinear双线性插值。不要看到一个库函数就用默认参数要搞清楚它的默认行为是什么不然数据上的偏差会很晦涩。重采样之后就可以分组统计了。假如你想统计每个县级行政区的夜间灯光总量、均值、最大值最直接的方案是用分区统计zonal statistics。这里我手动实现一份是希望你理解背后的逻辑后续哪怕换成超大数据集也心里有数。import geopandas as gpd import pandas as pd import numpy as np # 读取矢量边界保持与裁剪前一致的投影坐标系 polygons gpd.read_file(data/vector/study_area.shp) # 转换到与栅格数据相同的投影 polygons polygons.to_crs(resampled.rio.crs) # 用于提取像元中心坐标 from rasterio.features import geometry_mask # 这里简化先把栅格数据转成数组坐标 data resampled.values[0] # 单波段 xmin, ymin, xmax, ymax resampled.rio.bounds() width resampled.rio.width height resampled.rio.height # 像元大小 pixel_size_x (xmax - xmin) / width pixel_size_y (ymax - ymin) / height # 生成像元中心坐标网格 xs np.arange(xmin pixel_size_x/2, xmax, pixel_size_x) ys np.arange(ymax - pixel_size_y/2, ymin, -pixel_size_y) X, Y np.meshgrid(xs, ys)上面这段代码生成了像元中心坐标。接下来的思路是对每个矢量多边形用空间判断找出哪些像元落在里面。这里有很多现成库可以用比如rasterstats的zonal_stats但如果你数据量不大自己写也能练练手。更推荐的是直接装rasterstatsconda install -c conda-forge rasterstats -y用法非常简洁from rasterstats import zonal_stats stats zonal_stats( data/vector/study_area.shp, data/processed/02_重采样_500m.tif, stats[sum, mean, max, min, count], nodata0 ) # 转为 DataFrame 并和行政区名称合并 df_stats pd.DataFrame(stats) df_study gpd.read_file(data/vector/study_area.shp)[[name, geometry]] df_result pd.concat([df_study.drop(columnsgeometry), df_stats], axis1) print(df_result.head()) # 导出成 csv df_result.to_csv(data/processed/03_区域灯光统计.csv, indexFalse, encodingutf-8-sig)nodata0这里又一次出现是为了让统计函数跳过 NoData 像元。encodingutf-8-sig这个参数是我特意加的用 Excel 打开中文 CSV 时如果不用 utf-8-sig会直接乱码。老读者知道我有多少次倒在中文编码这个坑上这个参数可以说是国产数据从业者的必备细节。3.5 校验出图对比与数值抽查数据预处理做完不等于万事大吉。我始终觉得置信度比什么参数都重要所以一定要做校验。第一层校验是可视化把清洗、裁剪前后的影像画出来对比人工肉眼检查有没有明显错位或者空洞。fig, axes plt.subplots(1, 3, figsize(15, 5)) # 原始数据 with rasterio.open(tif_path) as src: original src.read(1, maskedTrue) axes[0].imshow(original, cmapcividis) axes[0].set_title(原始 NPP 夜间灯光) # 裁剪重采样后的数据 resampled_vis resampled.values[0] axes[1].imshow(resampled_vis, cmapcividis) axes[1].set_title(裁剪重采样后) # 热力图式展示 axes[2].imshow(resampled_vis, cmapmagma) axes[2].set_title(假彩色展示) plt.tight_layout() plt.show()如果影像边界与研究区矢量对不上比如整体偏移了几个像元那一定是投影或者裁剪步骤出了问题需要回头检查。第二层校验是数值抽查。从最终的统计表里挑出典型区域手动计算几个像元的均值看看和统计表是否一致。也可以在控制台里打印一些关键位置的像元值跟统计结果做交叉验证。对于一个汇总统计结果我还习惯检查三个环节总和是否为正、是否存在明显超界值、count 列的有效像元数是否合理。这些指标如果异常我就能在建模之前拦住错误而不是把脏数据送进下游流程。4. 踩过的坑与排查技巧实录最后这部分我把这些年实打实踩过的坑按主题列出来。这些坑并不高深但每一个都曾经浪费过我小半天时间。我把它们整理成一份“排查手册”希望你能绕开。4.1 坑一中文标签乱码图和表全废matplotlib 默认字体不支持中文这是老生常谈。但具体到 Notebook 里还有个隐蔽问题Notebook 输出的图在网页里能显示导出成图片却常常乱码。解决办法是在 Notebook 最开头就执行一次全局配置import matplotlib.pyplot as plt plt.rcParams[font.sans-serif] [SimHei, Microsoft YaHei, Arial Unicode MS] plt.rcParams[axes.unicode_minus] Falseaxes.unicode_minus False这行很多人会漏。它解决的是负号显示成方块的问题尤其在做数据预处理时经常遇到负值统计、差异图之类的图设置不设置差别非常大。另外如果上面这几个字体都没有也可以直接用plt.rcParams[font.family] sans-serif并指定你系统里已经安装的中文字体名字比如 Linux 下有Noto Sans CJK SC。4.2 坑二Notebook 跑着跑着内存爆了栅格数据看起来不大可一旦读到内存里尤其是多次读取、反复重采样之后很容易把 16G 内存占满。我有一次处理 30 米分辨率的影像原始文件才 800M裁剪重采样之后反而占了几 GB 内存Notebook 直接卡死。后来我总结出三条原则。第一条能用rasterio窗口读就按窗口读不要一张图整个载入。比如统计时可以用windowed reading分块处理。第二条每次跑完大计算手动清理一下大变量del resampled_vis然后调用gc.collect()不要指望内存自动释放。第三条中间结果及时落盘内存里只保留当前步骤必要的数据。尤其是你在 Notebook 里反复点击“重新运行”旧变量不会自动消失积少成多就会爆内存。建议每次修改代码前先Restart Kernel and Run All确保从干净的内核状态开始跑。4.3 坑三pandas 的 SettingWithCopyWarning处理表格数据时最常见的告警就是这个SettingWithCopyWarning。它的本质是你从一个 DataFrame 切片出子集然后对子集赋值但 pandas 无法判断你到底是改了原始 DataFrame 还是只改了一个副本。这个告警看着不致命但它经常意味着后面统计里的数据根本没有被真的改掉。我现在的习惯是凡是需要写回的子集一律用.copy()显式复制subset df[df[region] 华东].copy() subset[normalized] subset[night_light] / subset[night_light].max()不要为了省那点内存省略.copy()。数据预处理阶段安全性永远优先于性能。还有一种变体是链式索引df[df[a]1][b] 0这种写法我基本不用因为它不仅会触发告警而且赋值行为在不同 pandas 版本下都可能不同。4.4 坑四GDAL/Rasterio 环境冲突这个前面提过一些这里展开讲。如果你用 pip 安装 rasterio大概率会遇到底层 GDAL 的链接问题最常见的报错是ERROR 4: Unable to open ...或者OSError: Cant load GDAL library。根本原因是系统里有多个 GDAL 版本rasterio 加载了错误的那一个。我的经验是在 Windows 上用 conda 创建环境后rasterio、fiona、pyproj、geopandas 统一用 conda 安装不要混着 pip。如果你确实需要用 pip 装一些 conda 里没有的包装完再conda install --force-reinstall rasterio来修复底层库的链接关系。在 Linux 服务器上尽量避免自己从源码编译 GDAL直接conda install -c conda-forge gdal rasterio是最稳妥的。还有一种常见情况是rio.to_raster()在保存 tif 时提示Non-Georeferenced dataset。这种通常是因为内存中的数组已经丢掉了坐标系元数据或者你在操作时误把rio.write_crs()覆盖成 None。解决办法是在保存前重新声明resampled.rio.write_crs(resampled.coords[spatial_ref].attrs[crs_wkt], inplaceTrue) resampled.rio.to_raster(output.tif)4.5 坑五小细节但致命的类型问题这一条不算坑更像血泪经验总结。在数据预处理时pandas 的列经常读出来是 object 类型看着是数字却不能参与运算需要先pd.to_numeric转换。日期时间列也可能被读成字符串导致排序错乱。如果你是在做栅格和矢量的连接分析还要注意字段名是否匹配、空间索引是否建立。很多时候预处理代码本身没问题问题出在读进来的字段类型和预期不一致。我每次进入下一步操作前常用一行代码做检查print(df.dtypes)养成这类习惯之后你会发现自己写代码的速度并没有变快但是返工率会明显下降。4.6 常见问题速查表为了看着方便我把高频问题整理成一张速查表放在这里供你直接参照。问题现象常见原因首选排查路径Notebook 默认目录不对配置文件未生效或配置文件路径错误运行jupyter notebook --generate-config后检查配置路径中文出现方块或乱码matplotlib 字体配置缺失设置plt.rcParams[font.sans-serif]和axes.unicode_minus栅格读取报错无法打开GDAL 库冲突或路径错误确认路径存在用 conda 统一安装 gdal、rasterio统计结果偏差极大NoData 被当成有效值参与运算检查nodata与masked设置确保统计时跳过无效像元裁剪出空白影像矢量与栅格坐标系不一致输出前打印两边的 crs进行to_crsCSV 中文 Excel 打开乱码文件编码不是 UTF-8 with BOM导出时用encodingutf-8-sig内存溢出被内核杀大变量未释放、全图载入删除大变量并gc.collect()用分块读取或小范围裁剪5. 未完把自己变成流水线的一部分如果你能坚持读完前面这些内容说明你已经做好和“脏数据”长期共存的准备了。我个人的体会是数据预处理真正深奥的地方不是某个库的某个函数怎么写而是如何在每一步操作里保留可回溯的证据让你和任何一个后来接手的人都能说清楚“这份数据为什么被处理成了这样”。这也是我为什么把每次预处理的流程固定成步骤、固定成命名、固定成校验动作。技巧会过时库的 API 会更新但这种对过程透明性的坚持放在任何一个数据项目里都不会错。最后再分享一个小技巧每次开始一项新的预处理任务时我会在 Notebook 第一格写上项目名称、数据来源、创建日期、处理人和一句“处理目标”然后才碰任何代码。这些信息会随 Notebook 一起被保存、分享哪怕三个月后有人问起这个表格是怎么来的也能直接翻出答案一清二楚。很多时候数据工作者最重要的产出物不是模型精度的提升而是让处理过程经得起追问。希望这篇经验之谈对你有用。 Jupyter Notebook 数据预处理实战：从 NPP 夜间灯光 tif 到统计表
1. 为什么我坚持用 Jupyter Notebook 做数据预处理先交代一下背景。我日常工作中有一大半时间花在“把乱七八糟的数据弄成能用的数据”——这一步有个专业名字叫数据预处理。以前我习惯用 PyCharm 写一个完整的 .py 文件跑完看结果错了一步就得从头看日志改一处变量又要重新跑一遍全流程。直到我尝试把整个流程搬到 Jupyter Notebook 上才发现自己以前浪费了多少时间。如果在看这篇文章的你也是个经常和数据打交道的人我建议你认真看完因为数据预处理场景下交互式执行的优势是常态化的刚需而不是锦上添花。Jupyter Notebook 的核心价值可以概括为三个词分段执行、状态保留、可视化反馈。分段执行意味着我可以把“读数据—清洗—转换—统计—导出”拆成一个一个的单元格跑到哪里断了就在哪里修不用每次从头跑。状态保留则解决了“内存里的东西不丢”的问题你上一步读入的 DataFrame、计算好的掩膜、调整过的参数全部存在内核里后续单元格能直接引用。可视化反馈更是杀手锏一个matplotlib出图放在那里数据长什么样一目了然比对着 10 万行的表格看头部几行要直观得多。这套工作流不仅适合表格数据对栅格数据、遥感影像同样适用。我接下来会用一个 NPP 夜间灯光数据的例子完整走一遍数据预处理的各个环节把安装配置、清洗逻辑、裁剪投影、统计导出、踩坑排查都串起来。你可以把它当成一份直接抄的作业也可以在此基础上替换成你自己的数据。1.1 为什么交互式执行对预处理这么重要说白了数据预处理是个“试错”的过程。你拿到一份数据不知道它的缺失值长什么样、单位是什么、坐标系有没有偏离只能一个个条件去试。脚本语言里一次完整的python script.py执行通常要等几十秒甚至几分钟才能看到结果中间一个字段名写错又要重来。Notebook 则把整个过程变得像“聊天”一样写一个单元格、跑一个单元格立刻能看到输出。我现在的习惯是第一轮先建一个“侦察型 Notebook”专门用来摸数据的底比如 df.info()、df.describe()、df.head()然后是栅格数据的 crs、分辨率、nodata、波段数。这些探索性的判断通常会产生很多临时变量如果在脚本里最后要么删掉要么注释掉非常繁琐在 Notebook 里就无所谓留着也不碍事。第二轮再根据侦察结果写正式的清洗流程。这种“先侦察、后动刀”的思路能减少大量无效代码。还有一个实验室场景的好处就是可以保留实验的“中间状态”。我在处理夜间灯光数据的时候重投影、裁剪这类栅格操作非常耗时经常跑到一半我要去确认某个参数对不对。在 Notebook 里我可以随时暂停改掉参数然后只重跑当前单元格和后续部分前面的数据读入、预处理都不用重新执行。1.2 环境搭建的两种姿势与选择逻辑聊环境总有人问到底是装 Anaconda 还是 Miniconda还是直接用 pip 一个个装。我的看法是如果你不是已经有 Python 基础的老手直接上 Anaconda 最省心。因为 Anaconda 自带 Python 解释器、Jupyter Notebook 和一大堆常用库装完就能用而 Miniconda 只是最小化环境管理器还得手动创建环境和安装包。不要被很多教程里“Anaconda 太笨重”的说法吓到这年头磁盘空间没这么金贵把时间花在安装配置上才是真的不划算。如果一定要用命令行我的推荐是这么一条# 创建独立环境尽量指定 Python 版本避免日后依赖冲突 conda create -n geodata python3.9 -y conda activate geodata # 安装核心数据分析库 conda install -y jupyter pandas numpy matplotlib # 处理栅格、矢量数据的常用库 conda install -y rasterio geopandas rioxarray这里有个细节我明确要求python3.9而不是直接用默认的 Python 3.12 之类。原因很简单rasterio、GDAL 这类底层库的编译依赖非常重虽然官方维护者一直在跟新版本但偶尔还是会遇到某个轮子没跟上最新版 Python 的情况。用 3.9 或者 3.10在 conda-forge 上几乎能找到所有你需要的预编译包少很多麻烦。如果你是 pip 用户也有一条对应的路径pip install notebook pandas numpy matplotlib pip install rasterio geopandas rioxarray不过我个人的经验是GDAL 相关的生态库用 conda 安装的稳定性明显更好。因为 conda 不仅装 Python 包还会统一管理 GDAL、PROJ、GEOS 这些底层 C 库的版本。rasterio 这类库对底层库版本极度敏感pip 安装经常会遇到一个包依赖新版 GDAL、另一个包锁定旧版 GDAL 的冲突。conda 在这方面像个大管家会把所有底层依赖一起协调好。1.3 先装顺手的三件套目录、默认路径和快捷键很多人用 Jupyter Notebook都是打开浏览器新建一个 ipynb然后开始写。但用了一段时间你就会发现不解决三个问题Notebook 用起来会很难受一是文件多了之后目录乱成一团二是不知道文件默认保存到哪里来回切换烦人三是重复操作全靠鼠标效率很低。第一件事是给 Notebook 添加目录。老版本的 Jupyter Notebook比如 6.x可以通过插件实现命令是conda install -c conda-forge jupyter_contrib_nbextensions -y jupyter nbextension enable --py jupyter_contrib_nbextensions --sys-prefix装完重启 Jupyter打开任意一个 ipynb上方菜单栏如果有 Nbextensions 标签进去勾选 Table of Contents (2)页面上就会出现侧边目录。但注意如果你用的是 JupyterLab现在很多新版本的默认界面就是 Lab那么其实不用装任何插件它自带 Table of Contents 面板点左侧那个带层叠图标的面板就能看到当前 notebook 的目录结构。所以如果你在旧版 Notebook 里折腾了半天插件不生效我建议直接切换到 JupyterLab 试试。第二件事修改默认保存路径。默认情况下Jupyter Notebook 会把文件保存在启动时所在目录如果你每次都是从不同的地方启动文件就会散落各处。我自己的做法是把所有 ipynb 统一放在D:/jupyter_workspace下面操作方式分两步# 先生成配置文件 jupyter notebook --generate-config然后用文本编辑器打开生成的文件一般在用户目录下的.jupyter/jupyter_notebook_config.py找到这一行# c.NotebookApp.notebook_dir 改成c.NotebookApp.notebook_dir D:/jupyter_workspace保存之后再启动 Jupyter Notebook你会发现默认工作目录已经变了。这个配置也适用于 JupyterLab配置文件里改成c.ServerApp.notebook_dir也可以看你使用的版本。还有一个小技巧如果你只是想临时指定某个目录不需要改配置可以用jupyter notebook --notebook-dirD:/jupyter_workspace这种方式。第三件事是几个高频快捷键。在 Notebook 里按一次Esc进入命令模式然后B在当前单元格下方新建一个单元格A在当前单元格上方新建一个单元格D D连续按两次 D 删除当前单元格M把当前单元格切换为 MarkdownY把当前单元格切换为代码Shift Enter运行当前单元格并跳到下一个这组快捷键里我最常用的是Shift Enter和M。写数据预处理的文档时我会频繁插入 Markdown 分隔说明用快捷键切换模式比鼠标点击快得多。我的个人心得是数据预处理 Notebook 里Markdown 单元格的注释几乎和代码一样重要你一周后再回头看还能回忆起每一步在处理什么。2. 数据预处理不是写代码是先画流程图很多人处理数据的习惯是拿到数据直接开始写 pandas 操作写一步看一步最后数据处理完了整个流程却说不清楚。这种做法在小数据集上还行一旦数据量上来、步骤变多就很容易出现“处理完发现某个条件判断错了、前面所有清洗都白做”的悲剧。我现在的习惯是动手前先在 Notebook 顶部用 Markdown 画一个流程草图哪怕是几个列表也够用。这一步表面上看耽误了五分钟实际上能帮你节省至少半小时的返工时间。数据预处理说到底是一条流水线原始数据进来经过清洗、变换、合并、采样、校验最终变成可以进入建模或分析的规范数据集。流水线的每一道工序都要想清楚它要解决什么问题。2.1 一次完整预处理的六个步骤如果把数据预处理分成标准步骤大致是这六段数据加载与侦察读入数据了解字段类型、缺失值、分布、坐标系、分辨率等基本属性。数据清洗处理缺失值、重复值、异常值、非法值让数据在“物理层面”变得规范。数据变换重命名、类型转换、归一化、维度调整、重投影、重采样等让数据在“结构层面”变得可用。数据整合多张表拼接、栅格与矢量叠加、时间序列对齐等。数据统计与样本梳理生成描述性统计抽样核对数据合理性。数据导出与文档输出规范文件记录处理参数、版本号和步骤索引。这个框架不只是我的习惯几乎所有成熟的数据团队都有类似的分层。比如说你导出的最终数据出了问题你要能够回溯到它是从第几步开始坏的。如果中间产物都留着排查就很简单如果全都混在一个变量里覆盖来覆盖去就只能从头再来。2.2 把“中间产物规范”定下来我强烈建议在开始预处理之前先定义一个“产物规范”。这里指的是一套固定的命名、格式和组织规则比如原始数据单独放在raw/目录永远不修改修改后的数据放在processed/。中间产物文件名带上步骤号比如01_读入检查.pkl、02_清洗后.csv、03_裁剪.tif、04_重采样.tif。统一用YYYYMMDD作为日期后缀方便按版本回滚。统计结果和清洗后的数据分开存放不要混在一个文件里。这套规范看起来很简单实际操作中却能救你很多次。举个例子我在处理夜间灯光数据时因为中间要经过裁剪、重投影两次栅格运算每步生成的 tif 都很大运行时间也不短。如果我把每一步结果保存成一个中间 tif就算后面的操作出了问题我也只需要从对应的中间结果重新开始而不是把前面的运行又跑一遍。数据预处理有个铁律凡是跑得慢或不好恢复的步骤都值得做中间落盘。2.3 为什么选 NPP 夜间灯光数据做演示案例我做技术选型的时候喜欢挑一个“能覆盖绝大多数问题”的案例数据。NPP 夜间灯光数据恰好就是这么个典型。它是一种栅格影像数据常被用于城市扩展分析、区域社会发展研究等场景优点是数据相对容易获取、覆盖范围广、结构清晰缺点也很明显坐标系需要严格处理、原始数据有大量边界异常值、空间分辨率多种多样、文件体积不小。它的技术难点集中在四个方向一是栅格影像的读写不是寻常的 txt、csv需要 rasterio 这类专业库二是坐标系投影需要特殊处理用错投影直接导致空间偏移三是负值、极值、NoData 这类“脏数据”在影像里比表格里更隐蔽必须通过可视化才能发现四是栅格的统计方式与表格不同需要掩膜、分区统计等专业操作。任何一个玩过遥感或地理数据分析的人都会被这四个问题打磨过一遍。如果你只想快速体验流程可以用任何一张 GeoTIFF 栅格数据替代 NPP 夜间灯光数据步骤完全通用。但如果你想把整个案例完整复现那就跟着我的流程走你会发现这远比单纯处理一个 csv 学到的东西多得多。3. 完整实操NPP 夜间灯光数据预处理从 tif 到统计表下面进入正题。我用一台普通 Windows 笔记本环境是 conda 创建的geodata核心库是rasterio、geopandas、rioxarray。整个案例的目标是把原始 NPP 夜间灯光 GeoTIFF 文件读取、清洗、裁剪、重投影最终输出一个研究区内的像元统计表并生成一张可读的成果图。这里要做一个说明NPP-VIIRS 夜间灯光数据是 NASA/NOAA 合成的月度平均夜间灯光强度产品单位一般为nW/cm2/sr文件格式是 GeoTIFF数据类型有浮点型和整型两种版本。它每个像元的值代表该位置的平均夜间灯光辐亮度值越大通常说明该区域夜间灯光越强。它不像普通 CSV 那样可以直接在 Excel 里打开检查必须借助专业工具。3.1 侦察数据读进来先把“底细”摸清我拿到一份栅格数据从来不会急着写统计代码而是先做一个“读入检查”。这段代码信息量非常大import rasterio import numpy as np import matplotlib.pyplot as plt from rasterio.plot import show # 数据路径自行替换成你自己的文件 tif_path data/raw/VNL_v2_npp_202203_avg_vcmslcfg.tif with rasterio.open(tif_path) as src: print(CRS:, src.crs) print(尺寸(宽,高):, src.width, src.height) print(分辨率:, src.res) print(波段数:, src.count) print(NoData:, src.nodata) print(数据范围:, src.bounds) # 读取第一个波段 band1 src.read(1) print(数据类型:, band1.dtype) print(像元值范围:, np.nanmin(band1), np.nanmax(band1)) # 快速看一眼直方图 plt.figure(figsize(10, 4)) plt.hist(band1[band1 0].ravel(), bins50, range(0, 200)) plt.title(夜间灯光像元值直方图(仅正值)) plt.show()这段代码的输出会告诉我几件事坐标系是什么很多时候 NPP 数据用的是 WGS 84 经纬度坐标分辨率是多少常见的是 15 弧秒大概相当于赤道附近 500 米数据类型的浮点还是整型NoData 用的什么标记。这些信息直接决定了后面的处理参数。很多人在这一步就开始犯错了。有些人直接src.read(1)之后拿着数据就统计结果把 NoData 当成了真实 0 值或者把坐标投影忽略掉直接把经纬度当成平面坐标来算面积和距离。我建议把这一步的打印结果保存下来作为后续处理文档的“数据档案”。3.2 清洗处理负值、NoData 与极值NPP 夜间灯光数据的“脏”主要体现在负值和极端值上。月度合成的 VNL 数据虽然经过了初步处理但依然可能存在少量负值像元这些负值并不是真实灯光而多半是定标、几何校正过程中的噪声产物。我的清洗原则是负值一律先置为零或 NoData不要直接当有效值参与统计。import numpy as np # 重新读取 with rasterio.open(tif_path) as src: band src.read(1) profile src.profile # 查看负值像元数 neg_mask band 0 print(负值像元个数:, np.sum(neg_mask)) print(负值占比: {:.6f}%.format(np.sum(neg_mask) / band.size * 100)) # 清洗负值置为 NoData 或 0 # 这里选择置为0相对于负值更有利于后续统计 band[neg_mask] 0 # 防止极值干扰做一个上限截断 extreme_mask band 200 print(超过200的像元个数:, np.sum(extreme_mask)) # 如果确认是异常值可以截断否则保留 # band[extreme_mask] 200这个上限截断需要谨慎使用。夜间灯光数据在某些高亮度区域比如大型机场、商业中心确实会有很高的值直接截断反而会损失信息。我这里的extreme_mask只是先看一眼数量并不代表一定要处理。实际操作中我会结合研究区域的实际情况判断如果是做全局统计分析且异常像元占比极小我会选择保留如果异常像元挤占直方图导致可视化严重失真我会考虑对数变换或截断。处理异常值没有统一答案核心原则是每一次清洗操作都要记录下来操作原因和参数。NoData 的处理相对标准一些。如果用rioxarray可以这样import rioxarray raster rioxarray.open_rasterio(tif_path, maskedTrue) print(NoData 处理后的数据维度:, raster.shape) # maskedTrue 会把 NoData 位置变成 NaN print(NaN 像元数:, int(raster.isnull().sum()))采用maskedTrue的好处是后续的统计函数可以用nansum、nanmean这类函数直接跳过无效像元不用自己手写掩膜。这一点在处理有大量海域或云覆盖区域的影像时特别重要。3.3 按研究区裁剪并统一坐标系接下来是数据预处理的大头裁剪。也就是只保留你关心区域的像元。如果研究区是一个行政区或一个自定义矩形边界裁剪有两种常见方式按矢量边界裁剪和按经纬度范围裁剪。我以按矢量边界裁剪为例这也是遥感领域最常用的方式。import geopandas as gpd import rasterio from rasterio.mask import mask as rio_mask # 读取研究区矢量边界这里用一份简化面数据 study_area gpd.read_file(data/vector/study_area.shp) # 确保矢量与栅格在同一坐标系 # 如果矢量是 WGS84 经纬度而栅格也是 WGS84则跳过 reproject # 如果不是请用 study_area study_area.to_crs(栅格crs) crs_raster rasterio.open(tif_path).crs study_area study_area.to_crs(crs_raster) # 按边界裁剪 with rasterio.open(tif_path) as src: out_image, out_transform rio_mask(src, study_area.geometry, cropTrue, nodata0) out_meta src.meta.copy() print(裁剪后尺寸:, out_image.shape) # 更新元数据 out_meta.update({ driver: GTiff, height: out_image.shape[1], width: out_image.shape[2], transform: out_transform, nodata: 0 }) # 保存中间结果 with rasterio.open(data/processed/01_裁剪_研究区.tif, w, **out_meta) as dest: dest.write(out_image)这段逻辑里有几个容易踩的坑。第一裁剪之前必须检查矢量与栅格的坐标系是否一致不一致要to_crs转换。第二rio_mask的cropTrue表示对边界进行最小矩形裁剪这样输出影像会比原始影像小很多既能提升后续处理速度也能避免不必要的内存占用。第三nodata0这个参数要和前面清洗策略保持一致如果前面把负值置为 0这里裁掉边界外的部分也用 0 填充整条流水线才一致。如果你只需要一个矩形范围也可以用经纬度范围裁剪代码要简单很多# 用经纬度范围裁剪 min_lon, min_lat 110.0, 30.0 max_lon, max_lat 112.0, 32.0 if raster.rio.crs is None or raster.rio.crs.to_epsg() ! 4326: raster raster.rio.reproject(EPSG:4326) clipped raster.rio.clip_box(minxmin_lon, minymin_lat, maxxmax_lon, maxymax_lat)我这里用的rioxarray它封装了 rasterio很多操作一行就能搞定代码可读性也更好。你如果习惯纯 rasterio那就按第一个示例来两个思路是一样的。3.4 重采样与统计从栅格到表格先同步一个概念重采样就是把不同分辨率的栅格统一到一个目标分辨率。为什么要做这步因为很多时候你要拿夜间灯光数据和别的数据做叠加分析比如人口格网、土地利用影像它们的分辨率可能是 100 米、250 米、1000 米跟 NPP 的 500 米左右不一样。不统一分辨率像元之间无法对齐统计就无从谈起。import rioxarray # 重采样到目标分辨率例如500米 target_resolution 500 clipped rioxarray.open_rasterio(data/processed/01_裁剪_研究区.tif, maskedTrue) # 计算目标网格尺寸 # 先用一个简单的估算方法得到新分辨率下的行列数 width int((clipped.rio.bounds()[2] - clipped.rio.bounds()[0]) / target_resolution) height int((clipped.rio.bounds()[3] - clipped.rio.bounds()[1]) / target_resolution) # 重采样使用平均值方法对夜间灯光这种连续型变量更合适 resampled clipped.rio.reproject( clipped.rio.crs, resolutiontarget_resolution, resamplingrioxarray.enums.Resampling.average ) print(重采样后 shape:, resampled.shape) # 保存重采样结果 resampled.rio.to_raster(data/processed/02_重采样_500m.tif)关于采样算法我简单说一下选择逻辑对连续型数值变量比如夜间灯光辐亮度、气温、降水优先用average平均值重采样因为它不会丢失总辐射信息对分类变量比如土地利用类型优先用nearest最邻近或mode众数因为分类值不能被平均出“疑似新类别”如果是要精确配准、不希望改变像元值的范围可以用bilinear双线性插值。不要看到一个库函数就用默认参数要搞清楚它的默认行为是什么不然数据上的偏差会很晦涩。重采样之后就可以分组统计了。假如你想统计每个县级行政区的夜间灯光总量、均值、最大值最直接的方案是用分区统计zonal statistics。这里我手动实现一份是希望你理解背后的逻辑后续哪怕换成超大数据集也心里有数。import geopandas as gpd import pandas as pd import numpy as np # 读取矢量边界保持与裁剪前一致的投影坐标系 polygons gpd.read_file(data/vector/study_area.shp) # 转换到与栅格数据相同的投影 polygons polygons.to_crs(resampled.rio.crs) # 用于提取像元中心坐标 from rasterio.features import geometry_mask # 这里简化先把栅格数据转成数组坐标 data resampled.values[0] # 单波段 xmin, ymin, xmax, ymax resampled.rio.bounds() width resampled.rio.width height resampled.rio.height # 像元大小 pixel_size_x (xmax - xmin) / width pixel_size_y (ymax - ymin) / height # 生成像元中心坐标网格 xs np.arange(xmin pixel_size_x/2, xmax, pixel_size_x) ys np.arange(ymax - pixel_size_y/2, ymin, -pixel_size_y) X, Y np.meshgrid(xs, ys)上面这段代码生成了像元中心坐标。接下来的思路是对每个矢量多边形用空间判断找出哪些像元落在里面。这里有很多现成库可以用比如rasterstats的zonal_stats但如果你数据量不大自己写也能练练手。更推荐的是直接装rasterstatsconda install -c conda-forge rasterstats -y用法非常简洁from rasterstats import zonal_stats stats zonal_stats( data/vector/study_area.shp, data/processed/02_重采样_500m.tif, stats[sum, mean, max, min, count], nodata0 ) # 转为 DataFrame 并和行政区名称合并 df_stats pd.DataFrame(stats) df_study gpd.read_file(data/vector/study_area.shp)[[name, geometry]] df_result pd.concat([df_study.drop(columnsgeometry), df_stats], axis1) print(df_result.head()) # 导出成 csv df_result.to_csv(data/processed/03_区域灯光统计.csv, indexFalse, encodingutf-8-sig)nodata0这里又一次出现是为了让统计函数跳过 NoData 像元。encodingutf-8-sig这个参数是我特意加的用 Excel 打开中文 CSV 时如果不用 utf-8-sig会直接乱码。老读者知道我有多少次倒在中文编码这个坑上这个参数可以说是国产数据从业者的必备细节。3.5 校验出图对比与数值抽查数据预处理做完不等于万事大吉。我始终觉得置信度比什么参数都重要所以一定要做校验。第一层校验是可视化把清洗、裁剪前后的影像画出来对比人工肉眼检查有没有明显错位或者空洞。fig, axes plt.subplots(1, 3, figsize(15, 5)) # 原始数据 with rasterio.open(tif_path) as src: original src.read(1, maskedTrue) axes[0].imshow(original, cmapcividis) axes[0].set_title(原始 NPP 夜间灯光) # 裁剪重采样后的数据 resampled_vis resampled.values[0] axes[1].imshow(resampled_vis, cmapcividis) axes[1].set_title(裁剪重采样后) # 热力图式展示 axes[2].imshow(resampled_vis, cmapmagma) axes[2].set_title(假彩色展示) plt.tight_layout() plt.show()如果影像边界与研究区矢量对不上比如整体偏移了几个像元那一定是投影或者裁剪步骤出了问题需要回头检查。第二层校验是数值抽查。从最终的统计表里挑出典型区域手动计算几个像元的均值看看和统计表是否一致。也可以在控制台里打印一些关键位置的像元值跟统计结果做交叉验证。对于一个汇总统计结果我还习惯检查三个环节总和是否为正、是否存在明显超界值、count 列的有效像元数是否合理。这些指标如果异常我就能在建模之前拦住错误而不是把脏数据送进下游流程。4. 踩过的坑与排查技巧实录最后这部分我把这些年实打实踩过的坑按主题列出来。这些坑并不高深但每一个都曾经浪费过我小半天时间。我把它们整理成一份“排查手册”希望你能绕开。4.1 坑一中文标签乱码图和表全废matplotlib 默认字体不支持中文这是老生常谈。但具体到 Notebook 里还有个隐蔽问题Notebook 输出的图在网页里能显示导出成图片却常常乱码。解决办法是在 Notebook 最开头就执行一次全局配置import matplotlib.pyplot as plt plt.rcParams[font.sans-serif] [SimHei, Microsoft YaHei, Arial Unicode MS] plt.rcParams[axes.unicode_minus] Falseaxes.unicode_minus False这行很多人会漏。它解决的是负号显示成方块的问题尤其在做数据预处理时经常遇到负值统计、差异图之类的图设置不设置差别非常大。另外如果上面这几个字体都没有也可以直接用plt.rcParams[font.family] sans-serif并指定你系统里已经安装的中文字体名字比如 Linux 下有Noto Sans CJK SC。4.2 坑二Notebook 跑着跑着内存爆了栅格数据看起来不大可一旦读到内存里尤其是多次读取、反复重采样之后很容易把 16G 内存占满。我有一次处理 30 米分辨率的影像原始文件才 800M裁剪重采样之后反而占了几 GB 内存Notebook 直接卡死。后来我总结出三条原则。第一条能用rasterio窗口读就按窗口读不要一张图整个载入。比如统计时可以用windowed reading分块处理。第二条每次跑完大计算手动清理一下大变量del resampled_vis然后调用gc.collect()不要指望内存自动释放。第三条中间结果及时落盘内存里只保留当前步骤必要的数据。尤其是你在 Notebook 里反复点击“重新运行”旧变量不会自动消失积少成多就会爆内存。建议每次修改代码前先Restart Kernel and Run All确保从干净的内核状态开始跑。4.3 坑三pandas 的 SettingWithCopyWarning处理表格数据时最常见的告警就是这个SettingWithCopyWarning。它的本质是你从一个 DataFrame 切片出子集然后对子集赋值但 pandas 无法判断你到底是改了原始 DataFrame 还是只改了一个副本。这个告警看着不致命但它经常意味着后面统计里的数据根本没有被真的改掉。我现在的习惯是凡是需要写回的子集一律用.copy()显式复制subset df[df[region] 华东].copy() subset[normalized] subset[night_light] / subset[night_light].max()不要为了省那点内存省略.copy()。数据预处理阶段安全性永远优先于性能。还有一种变体是链式索引df[df[a]1][b] 0这种写法我基本不用因为它不仅会触发告警而且赋值行为在不同 pandas 版本下都可能不同。4.4 坑四GDAL/Rasterio 环境冲突这个前面提过一些这里展开讲。如果你用 pip 安装 rasterio大概率会遇到底层 GDAL 的链接问题最常见的报错是ERROR 4: Unable to open ...或者OSError: Cant load GDAL library。根本原因是系统里有多个 GDAL 版本rasterio 加载了错误的那一个。我的经验是在 Windows 上用 conda 创建环境后rasterio、fiona、pyproj、geopandas 统一用 conda 安装不要混着 pip。如果你确实需要用 pip 装一些 conda 里没有的包装完再conda install --force-reinstall rasterio来修复底层库的链接关系。在 Linux 服务器上尽量避免自己从源码编译 GDAL直接conda install -c conda-forge gdal rasterio是最稳妥的。还有一种常见情况是rio.to_raster()在保存 tif 时提示Non-Georeferenced dataset。这种通常是因为内存中的数组已经丢掉了坐标系元数据或者你在操作时误把rio.write_crs()覆盖成 None。解决办法是在保存前重新声明resampled.rio.write_crs(resampled.coords[spatial_ref].attrs[crs_wkt], inplaceTrue) resampled.rio.to_raster(output.tif)4.5 坑五小细节但致命的类型问题这一条不算坑更像血泪经验总结。在数据预处理时pandas 的列经常读出来是 object 类型看着是数字却不能参与运算需要先pd.to_numeric转换。日期时间列也可能被读成字符串导致排序错乱。如果你是在做栅格和矢量的连接分析还要注意字段名是否匹配、空间索引是否建立。很多时候预处理代码本身没问题问题出在读进来的字段类型和预期不一致。我每次进入下一步操作前常用一行代码做检查print(df.dtypes)养成这类习惯之后你会发现自己写代码的速度并没有变快但是返工率会明显下降。4.6 常见问题速查表为了看着方便我把高频问题整理成一张速查表放在这里供你直接参照。问题现象常见原因首选排查路径Notebook 默认目录不对配置文件未生效或配置文件路径错误运行jupyter notebook --generate-config后检查配置路径中文出现方块或乱码matplotlib 字体配置缺失设置plt.rcParams[font.sans-serif]和axes.unicode_minus栅格读取报错无法打开GDAL 库冲突或路径错误确认路径存在用 conda 统一安装 gdal、rasterio统计结果偏差极大NoData 被当成有效值参与运算检查nodata与masked设置确保统计时跳过无效像元裁剪出空白影像矢量与栅格坐标系不一致输出前打印两边的 crs进行to_crsCSV 中文 Excel 打开乱码文件编码不是 UTF-8 with BOM导出时用encodingutf-8-sig内存溢出被内核杀大变量未释放、全图载入删除大变量并gc.collect()用分块读取或小范围裁剪5. 未完把自己变成流水线的一部分如果你能坚持读完前面这些内容说明你已经做好和“脏数据”长期共存的准备了。我个人的体会是数据预处理真正深奥的地方不是某个库的某个函数怎么写而是如何在每一步操作里保留可回溯的证据让你和任何一个后来接手的人都能说清楚“这份数据为什么被处理成了这样”。这也是我为什么把每次预处理的流程固定成步骤、固定成命名、固定成校验动作。技巧会过时库的 API 会更新但这种对过程透明性的坚持放在任何一个数据项目里都不会错。最后再分享一个小技巧每次开始一项新的预处理任务时我会在 Notebook 第一格写上项目名称、数据来源、创建日期、处理人和一句“处理目标”然后才碰任何代码。这些信息会随 Notebook 一起被保存、分享哪怕三个月后有人问起这个表格是怎么来的也能直接翻出答案一清二楚。很多时候数据工作者最重要的产出物不是模型精度的提升而是让处理过程经得起追问。希望这篇经验之谈对你有用。