使用 terra package R 投影点
Projecting points with terra package R
我需要在 terra 包中投影 longitude/latitude 坐标，但我认为它无法正常工作，因为我试图使用此投影从栅格中提取数据，但数据不正确提取正确。
这是我的 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))
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/lat 是相同的）。
为什么你认为它与投影有关？无论哪种方式，它似乎都适合我。
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, basename(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[-] ~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)
这并不总是可行的，我需要在幕后做一些优化。
我需要在 terra 包中投影 longitude/latitude 坐标，但我认为它无法正常工作，因为我试图使用此投影从栅格中提取数据，但数据不正确提取正确。
这是我的 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))
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/lat 是相同的）。
为什么你认为它与投影有关？无论哪种方式，它似乎都适合我。
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, basename(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[-] ~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)
这并不总是可行的，我需要在幕后做一些优化。