0

我正在尝试从位于https://nomads.ncep.noaa.gov/pub/data/nccf/com/hrrr/prod/hrrr.20210511/conus/的 grib 文件中提取点数据。

这些文件由 170 或 173 个层(变量)组成。我想使用 terra:extract 提取点处的单元格值。但是,提取失败(层之间有些东西搞砸了)。这是我想为每个变量提取数据的 lon/lat:

latlon_df <- structure(list(Lon = c(-103.289, -96.6735, -96.9041, -96.76864, 
-102.4694, -96.6814, -97.7504, -99.6754, -96.4802, -103.0007, 
-96.8897, -101.8539, -103.9717, -101.253, -99.1134, -96.5849, 
-98.0301, -99.9537, -99.4601, -99.7122, -103.8278, -98.931, -102.1081, 
-101.7162, -100.115, -101.3448, -100.7805, -103.5606, -96.5302, 
-99.4156, -103.281, -100.0063, -97.9928, -100.7208, -98.5289, 
-96.762, -96.9218, -97.1024, -103.3793, -101.0841, -102.6745, 
-96.9188, -97.5154, -100.7435, -98.6938), Lat = c(45.5194, 44.3099, 
43.0526, 44.3252, 45.5183, 43.7316, 45.6796, 45.4406, 44.7154, 
44.0006, 43.7687, 43.9599, 43.4737, 44.9875, 45.0292, 44.0867, 
45.5735, 44.9895, 44.5256, 43.5938, 43.7343, 45.7163, 45.9189, 
43.1672, 45.6716, 45.9154, 45.7963, 44.6783, 44.5073, 43.7982, 
43.3784, 44.2912, 43.3841, 43.2002, 44.8579, 43.5048, 43.5033, 
45.1055, 44.4245, 45.4167, 44.5643, 44.304, 45.2932, 43.5601, 
43.7321)), class = "data.frame", row.names = c(NA, -45L))

library(terra)

latlons <- terra::vect(latlon_df,geom=c('Lon','Lat'),crs="+proj=longlat")
latlons <- project(latlons,"+proj=lcc +lat_0=38.5 +lon_0=262.5 +lat_1=38.5 +lat_2=38.5 +x_0=0 +y_0=0 
+R=6371229 +units=m +no_defs")

grib_data <- terra::rast(path to grib downloaded grib file)
extracted_data <- terra::extract(grib_data,latlons)

extract_data 是一团糟,没有任何意义。如果我一次做一层效果很好(但需要更长的时间)。

terra::extract(grib_data[[layer]],latlons)[,2]

我必须做什么才能让多层提取物起作用?

4

0 回答 0