Python处理NetCDF?Xarray怎么切片?

wowwwai GIS研习社 · 工具流程与项目排障

引言：很多 GIS 同学在做气象、海洋、遥感格网数据分析时，第一步就会遇到“Python处理NetCDF?Xarray怎么切片?”这个问题：NetCDF 文件能打开，但不知道如何按经纬度、时间、变量或空间范围取出自己需要的数据。本文用一个实际工作流说明 Xarray 读取 NetCDF、查看维度、按时间切片、按经纬度切片、导出结果的完整方法。

背景：为什么 GIS 场景经常需要 Python处理NetCDF

NetCDF 是气象、海洋、生态、水文和遥感格网数据中非常常见的格式。它和普通 Shapefile、GeoJSON 不一样，通常不是简单的二维要素表，而是一个多维数组文件。

一个典型 NetCDF 文件可能包含这些内容：

• 变量：例如 temperature、precipitation、wind_speed。
• 时间维度：例如每天、每小时、每月。
• 空间维度：通常是 lat、lon，也可能是 latitude、longitude、x、y。
• 高度或层级维度：例如 pressure_level、depth。
• 元数据：单位、坐标系、时间编码、缺失值等。

GIS 读者最常见的需求不是“读取整个 NetCDF”，而是只取一个区域、一个时间段、一个变量。也就是说，Python处理NetCDF的核心问题经常就是：Xarray怎么切片。

原理：Xarray切片本质上是在多维坐标上筛选数据

Xarray 可以把 NetCDF 读成带标签的多维数组。这里的“带标签”很关键：你不必只用数组下标 [0, 10, 20] 去取值，而可以直接按 time、lat、lon 这样的维度名切片。

常见对象有两个：

• Dataset：类似一个文件或数据集，里面可以有多个变量。
• DataArray：单个变量的数据数组，例如某一个降水变量。

在 GIS 工作中，理解这三个方法最重要：

方法	用途	适用场景
.sel()	按坐标标签选择	按时间、经纬度、行政区外接矩形切片
.isel()	按整数索引选择	取第一个时次、前 10 行、指定数组位置
.where()	按条件筛选	筛选大于某阈值的像元，或构造掩膜

如果你的数据维度是 time, lat, lon，那么 Xarray切片通常就是：

data.sel(time="2024-01-01", lat=slice(40, 20), lon=slice(100, 120))

注意，纬度切片方向要看文件中的 lat 是升序还是降序。这是 Python处理NetCDF 时最容易踩坑的地方之一。

步骤：用 Xarray 读取 NetCDF 并完成基础切片

1. 安装必要 Python 库

建议在独立环境中安装，避免和 QGIS、ArcGIS Pro 自带 Python 环境混用。

pip install xarray netcdf4 h5netcdf pandas numpy rioxarray rasterio

如果你只做 NetCDF 读取和切片，xarray、netcdf4 通常够用；如果后续要导出 GeoTIFF，建议安装 rioxarray 和 rasterio。

2. 打开 NetCDF 文件

import xarray as xr

nc_path = "data/temperature.nc"
ds = xr.open_dataset(nc_path)

print(ds)

输出信息中重点看四类内容：

• Dimensions：维度名称和长度。
• Coordinates：坐标变量，例如 time、lat、lon。
• Data variables：真正的数据变量。
• Attributes：单位、说明、坐标系等元数据。

例如你可能看到：

Dimensions:  (time: 365, lat: 181, lon: 360)
Coordinates:
  * time     (time) datetime64[ns] 2024-01-01 ... 2024-12-30
  * lat      (lat) float32 90.0 89.0 88.0 ... -89.0 -90.0
  * lon      (lon) float32 0.0 1.0 2.0 ... 358.0 359.0
Data variables:
    temp     (time, lat, lon) float32 ...

3. 选择一个变量

如果文件里有多个变量，先取出目标变量。比如温度变量名是 temp：

temp = ds["temp"]
print(temp)

这一步会得到一个 DataArray。后续的时间切片、经纬度切片都可以在这个对象上完成。

4. Xarray按时间切片

选择某一天：

temp_day = temp.sel(time="2024-01-15")

选择一个时间段：

temp_jan = temp.sel(time=slice("2024-01-01", "2024-01-31"))

如果时间维度不是标准日期，而是数字编码，先检查时间坐标：

print(ds["time"])
print(ds["time"].attrs)

多数符合 CF Convention 的 NetCDF 文件，Xarray 会自动解码时间。如果没有自动解码，可以打开时指定：

ds = xr.open_dataset(nc_path, decode_times=True)

5. Xarray按经纬度切片

假设要裁剪中国附近区域，经度 73 到 135，纬度 18 到 54。先检查纬度是升序还是降序：

print(temp["lat"].values[:5])
print(temp["lat"].values[-5:])

如果纬度从大到小，例如 90, 89, 88...，纬度切片应写成：

china = temp.sel(
    lon=slice(73, 135),
    lat=slice(54, 18)
)

如果纬度从小到大，例如 -90, -89, -88...，纬度切片应写成：

china = temp.sel(
    lon=slice(73, 135),
    lat=slice(18, 54)
)

这是 Xarray切片 和 GIS 软件裁剪工具的一个重要差异：slice() 的起止顺序要和坐标数组的实际排序一致。

6. 同时按时间和空间切片

实际项目中更常见的是同时筛选时间段和空间范围：

china_jan = temp.sel(
    time=slice("2024-01-01", "2024-01-31"),
    lon=slice(73, 135),
    lat=slice(54, 18)
)

如果变量维度不是 time, lat, lon 的顺序，也没有关系。Xarray 使用维度名选择，不依赖维度排列顺序。

7. 用最近邻方式选择一个站点附近格点

如果你有一个站点坐标，例如北京附近 lon=116.4, lat=39.9，可以选择最近的格点：

point_series = temp.sel(
    lon=116.4,
    lat=39.9,
    method="nearest"
)

print(point_series)

如果只取某一天的站点值：

point_value = temp.sel(
    time="2024-01-15",
    lon=116.4,
    lat=39.9,
    method="nearest"
)

print(float(point_value.values))

这在做气象站点匹配、模型验证、点位时间序列提取时非常实用。

8. 按整数索引切片：isel 的用法

如果你只是想快速查看第一个时次的数据，可以用 .isel()：

first_time = temp.isel(time=0)

取前 10 个时次：

first_10 = temp.isel(time=slice(0, 10))

.isel() 适合调试和快速抽样，但在正式 GIS 分析中，按日期和经纬度的 .sel() 更直观，也更不容易出错。

9. 将切片结果转成 DataFrame

如果结果是点位时间序列，可以转成 Pandas 表格：

df = point_series.to_dataframe(name="temp").reset_index()
print(df.head())

df.to_csv("output/beijing_temp_series.csv", index=False, encoding="utf-8-sig")

如果是区域格网数据，也可以转表，但要注意数据量可能非常大：

df_grid = china_jan.to_dataframe(name="temp").reset_index()
df_grid.to_csv("output/china_jan_temp_grid.csv", index=False, encoding="utf-8-sig")

对于大范围、多时间、多格点数据，不建议直接转 CSV。CSV 会迅速膨胀，后续读取也很慢。

10. 将切片结果导出为新的 NetCDF

如果后续还要保留时间维度和空间维度，推荐导出为 NetCDF：

china_jan.to_netcdf("output/china_jan_temp.nc")

这样可以继续被 Python、Panoply、QGIS、ArcGIS Pro 等工具读取。

11. 将单时次结果导出为 GeoTIFF

如果要在 QGIS 或 ArcGIS Pro 中制图，单时次二维格网更适合导出 GeoTIFF。下面示例使用 rioxarray：

import rioxarray

one_day = temp.sel(time="2024-01-15")

one_day = one_day.rio.set_spatial_dims(x_dim="lon", y_dim="lat")
one_day = one_day.rio.write_crs("EPSG:4326")

one_day.rio.to_raster("output/temp_20240115.tif")

如果你的坐标维度叫 longitude、latitude，需要对应修改：

one_day = one_day.rio.set_spatial_dims(x_dim="longitude", y_dim="latitude")

常见坑：Python处理NetCDF时 Xarray切片为什么会失败

1. 纬度方向写反，结果为空

最常见问题是这样的代码返回空数组：

subset = temp.sel(lat=slice(18, 54))

原因可能是文件里的纬度是从北到南排列，也就是 90, 89, 88...。这种情况下必须写成：

subset = temp.sel(lat=slice(54, 18))

排查方法：

print(temp["lat"].values[0], temp["lat"].values[-1])

2. 经度范围是 0 到 360，不是 -180 到 180

有些全球 NetCDF 的经度是 0 到 360，而你的 GIS 数据可能是 -180 到 180。例如美国西部经度 -120 在 0 到 360 体系中应为 240。

可以先查看经度范围：

print(float(temp["lon"].min()), float(temp["lon"].max()))

如果需要把 0 到 360 转成 -180 到 180，可参考：

ds = ds.assign_coords(lon=(((ds.lon + 180) % 360) - 180))
ds = ds.sortby("lon")

然后再按常规经纬度范围切片。

3. 维度名称不是 lat 和 lon

很多 NetCDF 文件不使用 lat、lon，而使用 latitude、longitude、x、y。不要凭经验写代码，先看结构：

print(ds.dims)
print(ds.coords)
print(ds.data_vars)

如果维度名是 latitude 和 longitude，切片代码应改为：

subset = temp.sel(
    longitude=slice(73, 135),
    latitude=slice(54, 18)
)

4. 时间没有被正确解码

如果 time 显示为整数而不是日期，可能是时间编码没有被正确解析。先查看属性：

print(ds["time"].attrs)

常见属性类似：

units: days since 1900-01-01
calendar: gregorian

如果自动解码失败，可以尝试：

ds = xr.open_dataset(nc_path, decode_times=True, use_cftime=True)

use_cftime 对非标准日历数据更友好，例如某些气候模式数据。

5. 数据太大，直接读取导致内存不足

NetCDF 可能非常大。Xarray 默认有延迟读取能力，但某些操作会触发真正加载。对于大文件，建议配合 Dask 分块：

ds = xr.open_dataset(
    nc_path,
    chunks={"time": 10}
)

如果同时按空间和时间切片，尽量先切片，再计算统计值：

subset = ds["temp"].sel(
    time=slice("2024-01-01", "2024-01-31"),
    lon=slice(73, 135),
    lat=slice(54, 18)
)

mean_temp = subset.mean(dim="time")

6. NetCDF 没有明确 CRS，导出 GeoTIFF 后位置不对

很多 NetCDF 数据只有经纬度坐标，没有像 GeoTIFF 那样显式写入 CRS。导出前要确认坐标是否真的是 WGS84 经纬度。如果是，则可写入：

data = data.rio.set_spatial_dims(x_dim="lon", y_dim="lat")
data = data.rio.write_crs("EPSG:4326")

如果数据是投影坐标，例如 Lambert Conformal Conic、极地投影或模式网格，不能简单写成 EPSG:4326，需要先正确识别原始坐标系。

方法比较：Xarray、Rasterio、GDAL、QGIS 处理 NetCDF 怎么选

方法	优势	局限	适合场景
Xarray	多维数据友好，按 time、lat、lon 切片直观	需要理解 Dataset、DataArray、维度和坐标	气象、海洋、气候 NetCDF 分析
Rasterio	GeoTIFF 读写方便，和栅格 GIS 工作流贴近	对多维 NetCDF 的时间维度处理不如 Xarray 直观	二维栅格读写、导出制图数据
GDAL	格式支持广，命令行批处理能力强	NetCDF 子数据集选择语法对新手不友好	格式转换、批量导出、服务器脚本
QGIS	可视化方便，适合检查数据位置和图层效果	复杂多维切片和批处理不如 Python 灵活	快速预览、制图、结果检查
ArcGIS Pro	多维栅格工具完整，适合桌面 GIS 工作流	自动化和开源生态灵活性较弱	企业 GIS 制图、多维栅格管理

如果你的目标是“批量提取某区域某时间段的 NetCDF 数据”，优先选 Xarray。如果目标是“做一张地图检查结果”，可以用 Xarray 切片后导出 GeoTIFF，再放到 QGIS 或 ArcGIS Pro 中检查。

检查清单：Xarray怎么切片才不容易出错

• 先用 print(ds) 查看 NetCDF 的完整结构，不要直接猜变量名。
• 确认目标变量名，例如 temp、t2m、precip。
• 确认空间维度名是 lat/lon、latitude/longitude，还是 x/y。
• 检查纬度是升序还是降序，再决定 slice() 的起止顺序。
• 检查经度范围是 -180 到 180 还是 0 到 360。
• 按时间切片前，确认 time 是否已被解码为日期。
• 大文件先做时间和空间裁剪，再计算均值、最大值、最小值。
• 导出 GeoTIFF 前，用 rioxarray 设置空间维度和 CRS。
• 导出后用 QGIS 或 ArcGIS Pro 检查位置、范围、像元值和 NoData。
• 不要把大范围多时次 NetCDF 直接转 CSV，除非明确知道数据规模。

FAQ：Python处理NetCDF与Xarray切片常见问题

Q1：Xarray怎么切片一个指定经纬度范围？

使用 .sel() 和 slice()。例如经度 73 到 135，纬度 18 到 54：

subset = temp.sel(
    lon=slice(73, 135),
    lat=slice(54, 18)
)

如果纬度是升序，则改为 lat=slice(18, 54)。关键是先检查纬度排序。

Q2：Python处理NetCDF时，为什么按 lat=slice(18,54) 没有结果？

大概率是纬度坐标为降序排列。Xarray 的 slice() 要跟坐标实际方向一致。如果 lat 从 90 到 -90，应使用 lat=slice(54,18)。

Q3：Xarray怎么按时间段切片？

如果时间已经被解析为日期，可以直接这样写：

subset = temp.sel(time=slice("2024-01-01", "2024-01-31"))

如果时间显示为数字，需要检查 time 的 units 和 calendar 属性，并尝试使用 decode_times=True 或 use_cftime=True。

Q4：NetCDF经度是0到360，怎么按中国范围切片？

中国范围经度 73 到 135 本身就在 0 到 360 范围内，可以直接切片。但如果你的研究区跨越负经度，例如 -120 到 -80，需要转换成 240 到 280，或者先把经度坐标统一转换到 -180 到 180。

Q5：Xarray切片后怎么导入 QGIS？

如果是单个时间片的二维格网，建议导出 GeoTIFF：

one_day = temp.sel(time="2024-01-15")
one_day = one_day.rio.set_spatial_dims(x_dim="lon", y_dim="lat")
one_day = one_day.rio.write_crs("EPSG:4326")
one_day.rio.to_raster("output/temp_20240115.tif")

然后在 QGIS 中直接加载这个 GeoTIFF，检查空间位置和数值范围。

Q6：Xarray的 sel 和 isel 有什么区别？

sel 按坐标标签选择，例如日期、经纬度；isel 按整数索引选择，例如第 0 个时间片、第 10 行。GIS 分析中，正式裁剪通常用 sel，调试抽样常用 isel。

Q7：多个 NetCDF 文件能不能一起切片？

可以。多个按时间分文件的 NetCDF 可以使用 open_mfdataset 合并读取：

ds = xr.open_mfdataset(
    "data/temp_2024_*.nc",
    combine="by_coords",
    chunks={"time": 10}
)

subset = ds["temp"].sel(
    time=slice("2024-01-01", "2024-01-31"),
    lon=slice(73, 135),
    lat=slice(54, 18)
)

前提是多个文件的变量名、坐标维度和元数据比较一致。

结论：先看维度，再用 sel 精准切片

Python处理NetCDF并不难，难点在于 NetCDF 是多维数据，不能用普通二维矢量或栅格的思路直接处理。Xarray 的优势是可以按 time、lat、lon 这些坐标标签进行清晰切片。

实际工作中建议固定一个流程：先 open_dataset，再 print(ds) 看变量和维度，接着检查经纬度方向和时间编码，最后使用 .sel() 完成时间与空间切片。需要制图时，再把单时次结果导出为 GeoTIFF，并在 QGIS 或 ArcGIS Pro 中检查。

只要把纬度方向、经度范围、时间解码和 CRS 这几个关键点检查清楚，Xarray切片就可以成为 GIS 项目中处理 NetCDF 数据的稳定基础工具。