Linux+Codex遥感卫星：农作物长势监测与地块划分实操精选
Sentinel-2 SAFE 数据包解压与完整性校验
先将 Sentinel-2 L2A 级 SAFE 压缩文件解包，进入 IMG_DATA/R10m 子目录，确认 B04（红光）与 B08（近红外）两个波段文件存在——这是后续 NDVI 计算的根基，缺一不可。流程虽然基础，但漏掉 B08 会导致 NDVI 无法生成，务必逐项核对。
在 Linux 服务器上处理遥感影像，核心思路是避开 Windows 专有 GIS 软件，全程使用开源命令行工具链。你已持有 Sentinel-2 L2A 级 SAFE 数据包，系统内安装了 GDAL 3.8+、Python 3.10+、rasterio、numpy、scikit-image 和 ogr2ogr——未装 QGIS 或 ArcGIS 正合此方案，纯命令行环境即可完成从数据加载、NDVI 计算、长势分级到地块矢量化输出的全流程。
切换到压缩包所在目录，执行：unzip S2B_MSIL2A_20260528T031549_N0509_R075_T49QEE_20260528T054220.zip。
解包后进入 SENTINEL2B_……/GRANULE/L2A_T49QEE_A029872_20260528T031549/IMG_DATA/R10m/ 子目录，检查是否存在 B02.jp2（蓝）、B03.jp2（绿）、B04.jp2（红）、B08.jp2（近红外）四个 10 米分辨率 JP2 文件——缺少 B08.jp2 将无法计算 NDVI，必须重新下载。
运行 gdalinfo B04.jp2 | grep "Size|Projection"，验证四个波段的图像尺寸一致（如 10980×10980），投影为 EPSG:32649（UTM Zone 49N），否则后续配准会出偏差。这一步不可跳过——少做一次校验，后续处理可能全部作废。
波段重采样与堆叠处理
方法一：用 gdalwarp 逐个重采样，统一分辨率与地理范围。执行：gdalwarp -tr 10 10 -t_srs EPSG:32649 -r bilinear -co COMPRESS=DEFLATE B02.jp2 B02_10m.tif，然后对 B03、B04、B08 重复操作，替换对应波段名即可。
方法二：一步堆叠更高效。先用 gdalbuildvrt -separate stack.vrt B02_10m.tif B03_10m.tif B04_10m.tif B08_10m.tif 构建虚拟栅格，再用 gdal_translate -co COMPRESS=DEFLATE stack.vrt stacked.tif 导出为 GeoTIFF。
需要注意：若原始 B08.jp2 本身是 20 米分辨率，强制重采样到 10 米会引入插值噪声。此时建议先用 gdalwarp -tr 20 20 将所有波段统一到 20 米，再整体缩放——盲目统一到 10 米会导致 NDVI 空间误差放大。这一细节常被忽视，但后果严重。
命令行计算 NDVI 并生成长势分级栅格
第一步：使用 rasterio + numpy 脚本计算 NDVI。新建 calc_ndvi.py，写入以下内容：
import rasterio, numpy as np
with rasterio.open("stacked.tif") as src:
red = src.read(3).astype(np.float32)
nir = src.read(4).astype(np.float32)
ndvi = (nir - red) / (nir + red + 1e-8)
profile = src.profile.copy()
profile.update(dtype=rasterio.float32, count=1, compress="DEFLATE")
with rasterio.open("ndvi.tif", "w", **profile) as dst:
dst.write(ndvi.astype(np.float32), 1)
第二步：若想快速验证结果而不写脚本，直接使用 gdal_calc.py：gdal_calc.py -A stacked.tif --A_band=3 -B stacked.tif --B_band=4 --outfile=ndvi_quick.tif --calc="(B-A)/(B+A+0.000001)" --type=Float32。两种方法任选其一。
第三步：按农业标准进行分级——这里采用非线性映射。先执行 gdal_translate -of GTiff -ot Byte -scale -1 0.8 1 255 ndvi.tif ndvi_byte.tif，再用 gdal_edit.py -a_nodata 0 ndvi_byte.tif 设置无效值。最后运行 python classify_ndvi.py ndvi_byte.tif longshi_classified.tif（需自行编写脚本，将 0–255 像素值映射为 1=较差、2=一般、3=良、4=优）。
基于 NDVI 阈值提取地块边界
① 使用 gdal_rasterize 将分类后的长势图转为矢量面，只提取“优”级区域：gdal_rasterize -burn 1 -where "DN=4" -l ndvi_classified ndvi_classified.tif good_zone.gpkg。
② 对矢量面进行形态学简化，消除细碎斑块：ogr2ogr -f GPKG -dialect sqlite -sql "SELECT ST_Buffer(geometry, 0), ST_Area(geometry) as area FROM good_zone WHERE area > 10000" good_buffered.gpkg good_zone.gpkg。面积小于 1 公顷的碎斑直接过滤，结果更干净。
③ 与行政边界叠加裁剪，确保最终结果落在真实耕地范围内：ogr2ogr -f GPKG -clipsrc county_boundary.shp final_good_fields.gpkg good_buffered.gpkg。命令看似简单，但存在关键陷阱：如果 county_boundary.shp 的坐标系不是 EPSG:32649，必须先用 ogr2ogr -t_srs EPSG:32649 重投影，否则裁剪结果完全错位。许多人在这步翻车——坐标系不统一带来的偏移可达几百米。