提取两个重叠栅格的数据
Extracting data of two overlapping rasters
我正在尝试提取两组重叠栅格的数据（一个是一组 35 个栅格，全部来自同一源，第二个是高程栅格）以获得 data.frame 个值（所有栅格的每个像素的值的平均值）。
光栅堆栈的描述如下：
> stack_pacifico
class : RasterStack
dimensions : 997, 709, 706873, 35 (nrow, ncol, ncell, nlayers)
resolution : 0.008333333, 0.008333333 (x, y)
extent : -81.62083, -75.7125, 0.3458336, 8.654167 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
names : F101992.v//ts.avg_vis, F101993.v//ts.avg_vis, F101994.v//ts.avg_vis, F121994.v//ts.avg_vis, F121995.v//ts.avg_vis, F121996.v//ts.avg_vis, F121997.v//ts.avg_vis, F121998.v//ts.avg_vis, F121999.v//ts.avg_vis, F141997.v//ts.avg_vis, F141998.v//ts.avg_vis, F141999.v//ts.avg_vis, F142000.v//ts.avg_vis, F142001.v//ts.avg_vis, F142002.v//ts.avg_vis, ...
min values : 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
max values : 63, 63, 63, 63, 63, 63, 63, 63, 63, 63, 63, 63, 63, 63, 63, ...
对于高程栅格：
> elevation_pacifico
class : RasterLayer
dimensions : 997, 709, 706873 (nrow, ncol, ncell)
resolution : 0.008333333, 0.008333333 (x, y)
extent : -81.62083, -75.7125, 0.3458336, 8.654167 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
data source : in memory
names : COL_msk_alt
values : -16, 5164 (min, max)
这是我第一次使用栅格数据，我想按 1km2（或更小）的网格提取数据。我知道可以强制两个栅格的分辨率以适应该区域要求，而且两个尺寸都相等，因此每个栅格的像素数相同。
我的问题是，我是否可以仅合并所有栅格（堆栈中的栅格和高程栅格）并在确信所有像素重叠（或位于同一位置）的情况下提取数据？或者我是否必须创建主 SpatialGrid 或 SpatialPixel 对象，然后将栅格数据提取到这些对象？
提前致谢，
点这里link可以下载光栅栈的数据（如果要下载所有的栈，可以使用https://github.com/ivanhigueram/nightlights中的脚本）：
http://www.ngdc.noaa.gov/eog/data/web_data/v4composites/
海拔：
#Download country map and filter by pacific states
colombia_departments <- getData("GADM", download=T, country="CO", level=1)
pacific_littoral <- c(11, 13, 21, 30)
pacific_littoral_map <- colombia_departments[colombia_departments@data$ID_1 %in% pacific_littoral, ]
#Download elevation data and filter it for pacific states
elevation <- getData("alt", co="COL")
elevation_pacifico <- crop(elevation, pacific_littoral_map)
elevation_pacifico <- setExtent(elevation_pacifico, rasters_extent)
如果两个栅格对象的分辨率、范围和坐标系完全相同，则单元格将完美重叠。您可以通过查看坐标来确认这一点：
coordinates(stack_pacifico)
coordinates(elevation_pacifico)
# are they the same?
identical(coordinates(stack_pacifico), coordinates(elevation_pacifico))
您可以使用以下方法之一提取每个对象的所有单元格值：
as.data.frame(r)
values(r)
r[]
extract(r, seq_len(ncell(r)))
其中 r 是您的光栅对象。
（对于单个栅格层，这些并不都具有一致的行为 - as.data.frame(r) 确保结果是 data.frame，如果 r 是单个栅格，则结果将只有一个列图层；相比之下，如果与单个栅格图层一起使用，替代方案将 return 是一个简单的矢量。）
as.data.frame(stack_pacifico) 的行对应于与 as.data.frame(elevation_pacifico) (or, equivalently, the elements ofvalues(elevation_pacifico)` 的行相同坐标的单元格。
或者这样做：
s <- stack(elevation_pacifico, stack_pacifico)
d <- values(s)
我正在尝试提取两组重叠栅格的数据（一个是一组 35 个栅格，全部来自同一源，第二个是高程栅格）以获得 data.frame 个值（所有栅格的每个像素的值的平均值）。
光栅堆栈的描述如下：
> stack_pacifico
class : RasterStack
dimensions : 997, 709, 706873, 35 (nrow, ncol, ncell, nlayers)
resolution : 0.008333333, 0.008333333 (x, y)
extent : -81.62083, -75.7125, 0.3458336, 8.654167 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
names : F101992.v//ts.avg_vis, F101993.v//ts.avg_vis, F101994.v//ts.avg_vis, F121994.v//ts.avg_vis, F121995.v//ts.avg_vis, F121996.v//ts.avg_vis, F121997.v//ts.avg_vis, F121998.v//ts.avg_vis, F121999.v//ts.avg_vis, F141997.v//ts.avg_vis, F141998.v//ts.avg_vis, F141999.v//ts.avg_vis, F142000.v//ts.avg_vis, F142001.v//ts.avg_vis, F142002.v//ts.avg_vis, ...
min values : 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
max values : 63, 63, 63, 63, 63, 63, 63, 63, 63, 63, 63, 63, 63, 63, 63, ...
对于高程栅格：
> elevation_pacifico
class : RasterLayer
dimensions : 997, 709, 706873 (nrow, ncol, ncell)
resolution : 0.008333333, 0.008333333 (x, y)
extent : -81.62083, -75.7125, 0.3458336, 8.654167 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
data source : in memory
names : COL_msk_alt
values : -16, 5164 (min, max)
这是我第一次使用栅格数据，我想按 1km2（或更小）的网格提取数据。我知道可以强制两个栅格的分辨率以适应该区域要求，而且两个尺寸都相等，因此每个栅格的像素数相同。
我的问题是，我是否可以仅合并所有栅格（堆栈中的栅格和高程栅格）并在确信所有像素重叠（或位于同一位置）的情况下提取数据？或者我是否必须创建主 SpatialGrid 或 SpatialPixel 对象，然后将栅格数据提取到这些对象？
提前致谢，
点这里link可以下载光栅栈的数据（如果要下载所有的栈，可以使用https://github.com/ivanhigueram/nightlights中的脚本）：
http://www.ngdc.noaa.gov/eog/data/web_data/v4composites/
海拔：
#Download country map and filter by pacific states
colombia_departments <- getData("GADM", download=T, country="CO", level=1)
pacific_littoral <- c(11, 13, 21, 30)
pacific_littoral_map <- colombia_departments[colombia_departments@data$ID_1 %in% pacific_littoral, ]
#Download elevation data and filter it for pacific states
elevation <- getData("alt", co="COL")
elevation_pacifico <- crop(elevation, pacific_littoral_map)
elevation_pacifico <- setExtent(elevation_pacifico, rasters_extent)
如果两个栅格对象的分辨率、范围和坐标系完全相同，则单元格将完美重叠。您可以通过查看坐标来确认这一点：
coordinates(stack_pacifico)
coordinates(elevation_pacifico)
# are they the same?
identical(coordinates(stack_pacifico), coordinates(elevation_pacifico))
您可以使用以下方法之一提取每个对象的所有单元格值：
as.data.frame(r)
values(r)
r[]
extract(r, seq_len(ncell(r)))
其中 r 是您的光栅对象。
（对于单个栅格层，这些并不都具有一致的行为 - as.data.frame(r) 确保结果是 data.frame，如果 r 是单个栅格，则结果将只有一个列图层；相比之下，如果与单个栅格图层一起使用，替代方案将 return 是一个简单的矢量。）
as.data.frame(stack_pacifico) 的行对应于与 as.data.frame(elevation_pacifico) (or, equivalently, the elements ofvalues(elevation_pacifico)` 的行相同坐标的单元格。
或者这样做：
s <- stack(elevation_pacifico, stack_pacifico)
d <- values(s)