在点位置的多个时间步长上提取 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 文件中每个时间步长的变量值(例如表面温度)。
我有一个 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 文件中每个时间步长的变量值(例如表面温度)。