问题描述
我需要在 terra 包中投影经度/纬度坐标，但我不相信它可以正常工作，因为我试图从具有此投影的栅格中提取数据，但未正确提取数据。
这是我的经纬度和我用来尝试投影它们的代码。
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))
latlons <- terra::vect(latlon_df,geom=c('Lon','Lat'),crs="+proj=longlat")
lcc <- terra::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")
var_df <- terra::extract(grib_data,lcc)[,-1]
我使用的栅格数据 (grib_data) 来自这里（它太大了，我无法放在这里）。 https://nomads.ncep.noaa.gov/pub/data/nccf/com/hrrr/prod/hrrr.20210612/conus/hrrr.t00z.wrfsubhf00.grib2
我不确定我在这里做错了什么，因为我以前使用过这种方法，而且它似乎工作正常。任何帮助都会很棒。
编辑：我遇到的具体问题是我没有为每个 lon/lat 对获得任何不同的值。每个变量的值不同，但所有站点的值（不同的lon/lats是相同的）。
解决方法
为什么你认为它与投影有关？无论哪种方式，它似乎都适合我。
url <- "https://nomads.ncep.noaa.gov/pub/data/nccf/com/hrrr/prod/hrrr.20210612/conus/hrrr.t00z.wrfsubhf00.grib2"
if (!file.exist(basename(url))) download.file(url,basename(url),mode="wb")
url <- paste0(url,".idx")
if (!file.exist(basename(url))) download.file(url,mode="wb")
library(terra)
r <- rast("hrrr.t00z.wrfsubhf00.grib2")
r
#class : SpatRaster
#dimensions : 1059,1799,49 (nrow,ncol,nlyr)
#resolution : 3000,3000 (x,y)
#extent : -2699020,2697980,-1588806,1588194 (xmin,xmax,ymin,ymax)
#coord. ref. : +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
#source : hrrr.t00z.wrfsubhf00.grib2
#names : 0[-] ~here",0[-] ~tops",0[-] ~here",0[-] ~face",1000[~ound",...
您可以检查与栅格数据重叠的点
plot(r,1)
points(lcc)
然后提取。 grib 文件需要很长时间，但它似乎确实有效
e <- extract(r,lcc)
head(e[,c(1,6,9)])
# ID 0[-] SFC="Ground or water surface" 0[-] SFC="Ground or water surface".1
#1 1 85100 11.775471
#2 2 54400 11.087971
#3 3 79300 9.900471
#4 4 49200 10.712971
#5 5 70800 9.212971
#6 6 56600 11.400471
确保您拥有当前的 (CRAN) 版本，或者您可以像这样安装的开发版本：
install.packages('terra',repos='https://rspatial.r-universe.dev')
您可以通过从磁盘进行单次读取来大大加快速度（在本例中通过添加零）
e <- extract(r+0,lcc)
这并不总是可行的，我需要在幕后做一些优化。