在点位置的多个时间步长上提取 netCDF 变量的值

Extract the values of netCDF variables over many time steps at point locations

我有一个 netCDF 文件,其中包含几个气象变量(例如边界层高度、地表温度、地表压力),在一个区域的 10 年期间每天测量一次。因此,每个变量有超过 4,000 个时间步长。

我想使用 QGIS 或 Python 在整个时间序列的几个点位置提取每个气象变量的值。所需的输出是 (a) 由点数据和提取的所有时间步长的气象变量值组成的新向量层,或 (b) 类似的输出,但采用 CSV/Excel 格式。

QGIS 中的 “Add raster values to points” 工具执行了所需的操作,但我无法将其应用于我的 NC 文件,因为它是网格数据,而不是栅格数据.这里的任何帮助将不胜感激。我也很满意 Python 是否有一个 Python 包可以实现这一点。

这是我的 netCDF 文件的 ncdump 输出:

netcdf Beijing_4UTC_2000-10 {
dimensions:
        longitude = 36 ;
        latitude = 19 ;
        time = 4018 ;
variables:
        float longitude(longitude) ;
                longitude:units = "degrees_east" ;
                longitude:long_name = "longitude" ;
        float latitude(latitude) ;
                latitude:units = "degrees_north" ;
                latitude:long_name = "latitude" ;
        int time(time) ;
                time:units = "hours since 1900-01-01 00:00:00.0" ;
                time:long_name = "time" ;
                time:calendar = "gregorian" ;
        short u100(time, latitude, longitude) ;
                u100:scale_factor = 0.000622911066002879 ;
                u100:add_offset = 0.118876376345661 ;
                u100:_FillValue = -32767s ;
                u100:missing_value = -32767s ;
                u100:units = "m s**-1" ;
                u100:long_name = "100 metre U wind component" ;
        short v100(time, latitude, longitude) ;
                v100:scale_factor = 0.000694160650027236 ;
                v100:add_offset = -3.10686248788727 ;
                v100:_FillValue = -32767s ;
                v100:missing_value = -32767s ;
                v100:units = "m s**-1" ;
                v100:long_name = "100 metre V wind component" ;
        short u10(time, latitude, longitude) ;
                u10:scale_factor = 0.000472531631472288 ;
                u10:add_offset = -0.729659788764955 ;
                u10:_FillValue = -32767s ;
                u10:missing_value = -32767s ;
                u10:units = "m s**-1" ;
                u10:long_name = "10 metre U wind component" ;
        short v10(time, latitude, longitude) ;
                v10:scale_factor = 0.000513328716515632 ;
                v10:add_offset = -3.10863126000035 ;
                v10:_FillValue = -32767s ;
                v10:missing_value = -32767s ;
                v10:units = "m s**-1" ;
                v10:long_name = "10 metre V wind component" ;
        short d2m(time, latitude, longitude) ;
                d2m:scale_factor = 0.00102728301940017 ;
                d2m:add_offset = 268.55166378769 ;
                d2m:_FillValue = -32767s ;
                d2m:missing_value = -32767s ;
                d2m:units = "K" ;
                d2m:long_name = "2 metre dewpoint temperature" ;
        short t2m(time, latitude, longitude) ;
                t2m:scale_factor = 0.00110899471613419 ;
                t2m:add_offset = 279.266550605181 ;
                t2m:_FillValue = -32767s ;
                t2m:missing_value = -32767s ;
                t2m:units = "K" ;
                t2m:long_name = "2 metre temperature" ;
        short blh(time, latitude, longitude) ;
                blh:scale_factor = 0.0884957109728695 ;
                blh:add_offset = 2911.24127842015 ;
                blh:_FillValue = -32767s ;
                blh:missing_value = -32767s ;
                blh:units = "m" ;
                blh:long_name = "Boundary layer height" ;
        short sp(time, latitude, longitude) ;
                sp:scale_factor = 0.374649985503487 ;
                sp:add_offset = 92661.2189250073 ;
                sp:_FillValue = -32767s ;
                sp:missing_value = -32767s ;
                sp:units = "Pa" ;
                sp:long_name = "Surface pressure" ;
                sp:standard_name = "surface_air_pressure" ;
        short tp(time, latitude, longitude) ;
                tp:scale_factor = 3.21524733289788e-07 ;
                tp:add_offset = 0.0105350794109732 ;
                tp:_FillValue = -32767s ;
                tp:missing_value = -32767s ;
                tp:units = "m" ;
                tp:long_name = "Total precipitation" ;

// global attributes:
                :Conventions = "CF-1.6" ;
                :history = "2020-06-30 13:19:05 GMT by grib_to_netcdf-2.16.0: /opt/ecmwf/eccodes/bin/grib_to_netcdf -S param -o /cache/data3/adaptor.mars.internal-1593521146.6122832-30945-36-b3534028-e955-439d-b3fd-a960d29aa305.nc /cache/tmp/b3534028-e955-439d-b3fd-a960d29aa305-adaptor.mars.internal-1593521146.6128585-30945-7-tmp.grib" ;
}

使用 ArcGIS Pro 中的“制作 NetCDF Table 视图”工具设法解决了这个问题。我选择了所有感兴趣的变量,将行维度设置为 'Time' 变量,并使用 lat/long 坐标作为维度值对。此 returns a table 具有在指定 XY 位置采样的 netCDF 文件中每个时间步长的变量值(例如表面温度)。