使用 extract() 函数在点处提取栅格数据:如何防止 NA 的输出
Extracting Raster Data at Points using the extract() function: How To Prevent the Output of NA's
我有一个 data.frame 具有三个变量(ID、经度和纬度)的海豚观察。
我还有 59 个 AQUA/MODIS 包含海面温度 (SST) 的 netCDF 文件。我想提取海豚观测位置的 SST 栅格数据。
海豚数据
d <- read.csv("dolphins.csv")
str(d)
# 'data.frame': 650 obs. of 4 variables:
#$ ID : num 1 2 3 4 5 6 7 8 9 10 ...
#$ Date : Factor w/ 73 levels "1/24/17","1/24/18",..: 67 67 67 8 8 8 8 8 8 8 ...
#$ Latitude : num 42 42 42 42 42 ...
#$ Longitude: num 19.1 19.1 19.1 19.1 19.1 ...
SST 栅格
library(terra)
filenames = list.files('Ocean_ColorSST_2016_2021',pattern='*.nc',full.names=TRUE)
SSTs <- rast(filenames, "sst")
SSTs
#class : SpatRaster
#dimensions : 4320, 8640, 59 (nrow, ncol, nlyr)
#resolution : 0.04166667, 0.04166667 (x, y)
#extent : -180, 180, -90.00001, 90 (xmin, xmax, ymin, ymax)
#coord. ref. : lon/lat WGS 84
#sources : AQUA_MODIS.20160901_20160930.L3m.MO.SST.sst.4km.nc:sst
# AQUA_MODIS.20161001_20161031.L3m.MO.SST.sst.4km.nc:sst
# AQUA_MODIS.20161101_20161130.L3m.MO.SST.sst.4km.nc:sst
# ... and 56 more source(s)
#varnames : sst (Sea Surface Temperature)
# sst (Sea Surface Temperature)
# sst (Sea Surface Temperature)
# ...
#names : sst, sst, sst, sst, sst, sst, ...
#unit : degree_C, degree_C, degree_C, degree_C, degree_C, degree_C, ...
创建一个 SpatialPointsDataFrame
library(sp)
points_spdf_M <- d
coordinates(points_spdf_M) <- ~ Latitude + Longitude
crs(points_spdf_M) <- crs(SSTs)
points_spdf_M
#class : SpatialPoints
#features : 650
#extent : 41.5978, 42.67778, 18.24337, 19.99933 (xmin, xmax, ymin, ymax)
#crs : +proj=longlat +datum=WGS84 +no_defs
从海豚 ID returns NAs
中提取光栅数据
library(raster)
ncin_SST <- stack(SSTs)
Extract_SST_M <- extract(ncin_SST, points_spdf_M)
head(Extract_SST_M)
#sst.1 sst.2 sst.3 sst.4 sst.5 sst.6 sst.7 sst.8 sst.9 sst.10 sst.11 sst.12 sst.13 sst.14 sst.15 sst.16 sst.17 sst.18
# [1,] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
# [2,] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
# [3,] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
我试过的另一种方法
v <- vect(d, geom=c("Longitude", "Latitude"))
e <- extract(SSTs, v)
head(e)
输出 ?
你像这样创建一个 SpatialPointsDataFrame
coordinates(Final_M_Points) <- ~ Latitude + Longitude
它应该在哪里
coordinates(Final_M_Points) <- ~ Longitude + Latitude
您的工作流程得到简化 --- 仅使用 terra
library(terra)
filenames <- list.files("---")
SSTs <- rast(filenames, "sst")
d <- read.csv("file.csv")
v <- vect(d, geom=c("Longitude", "Latitude"))
e <- extract(SSTs, v)
最好用terra
,避免sp
和raster
。在这种情况下,您会收到经度和纬度颠倒的警告。
我有一个 data.frame 具有三个变量(ID、经度和纬度)的海豚观察。
我还有 59 个 AQUA/MODIS 包含海面温度 (SST) 的 netCDF 文件。我想提取海豚观测位置的 SST 栅格数据。
海豚数据
d <- read.csv("dolphins.csv")
str(d)
# 'data.frame': 650 obs. of 4 variables:
#$ ID : num 1 2 3 4 5 6 7 8 9 10 ...
#$ Date : Factor w/ 73 levels "1/24/17","1/24/18",..: 67 67 67 8 8 8 8 8 8 8 ...
#$ Latitude : num 42 42 42 42 42 ...
#$ Longitude: num 19.1 19.1 19.1 19.1 19.1 ...
SST 栅格
library(terra)
filenames = list.files('Ocean_ColorSST_2016_2021',pattern='*.nc',full.names=TRUE)
SSTs <- rast(filenames, "sst")
SSTs
#class : SpatRaster
#dimensions : 4320, 8640, 59 (nrow, ncol, nlyr)
#resolution : 0.04166667, 0.04166667 (x, y)
#extent : -180, 180, -90.00001, 90 (xmin, xmax, ymin, ymax)
#coord. ref. : lon/lat WGS 84
#sources : AQUA_MODIS.20160901_20160930.L3m.MO.SST.sst.4km.nc:sst
# AQUA_MODIS.20161001_20161031.L3m.MO.SST.sst.4km.nc:sst
# AQUA_MODIS.20161101_20161130.L3m.MO.SST.sst.4km.nc:sst
# ... and 56 more source(s)
#varnames : sst (Sea Surface Temperature)
# sst (Sea Surface Temperature)
# sst (Sea Surface Temperature)
# ...
#names : sst, sst, sst, sst, sst, sst, ...
#unit : degree_C, degree_C, degree_C, degree_C, degree_C, degree_C, ...
创建一个 SpatialPointsDataFrame
library(sp)
points_spdf_M <- d
coordinates(points_spdf_M) <- ~ Latitude + Longitude
crs(points_spdf_M) <- crs(SSTs)
points_spdf_M
#class : SpatialPoints
#features : 650
#extent : 41.5978, 42.67778, 18.24337, 19.99933 (xmin, xmax, ymin, ymax)
#crs : +proj=longlat +datum=WGS84 +no_defs
从海豚 ID returns NAs
中提取光栅数据library(raster)
ncin_SST <- stack(SSTs)
Extract_SST_M <- extract(ncin_SST, points_spdf_M)
head(Extract_SST_M)
#sst.1 sst.2 sst.3 sst.4 sst.5 sst.6 sst.7 sst.8 sst.9 sst.10 sst.11 sst.12 sst.13 sst.14 sst.15 sst.16 sst.17 sst.18
# [1,] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
# [2,] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
# [3,] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
我试过的另一种方法
v <- vect(d, geom=c("Longitude", "Latitude"))
e <- extract(SSTs, v)
head(e)
输出 ?
你像这样创建一个 SpatialPointsDataFrame
coordinates(Final_M_Points) <- ~ Latitude + Longitude
它应该在哪里
coordinates(Final_M_Points) <- ~ Longitude + Latitude
您的工作流程得到简化 --- 仅使用 terra
library(terra)
filenames <- list.files("---")
SSTs <- rast(filenames, "sst")
d <- read.csv("file.csv")
v <- vect(d, geom=c("Longitude", "Latitude"))
e <- extract(SSTs, v)
最好用terra
,避免sp
和raster
。在这种情况下,您会收到经度和纬度颠倒的警告。