【问题标题】:Extracting data of two overlapping rasters提取两个重叠栅格的数据
【发布时间】:2016-01-08 19:56:14
【问题描述】:

我正在尝试提取两组重叠栅格的数据(一个是 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 对象,然后将栅格数据提取到这些对象?

提前致谢,

点击此链接可以下载光栅栈中的数据(如果要下载所有的栈,可以使用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)

【问题讨论】:

  • 如果两个堆栈的分辨率、范围和坐标系相同,则单元格将完美重叠。 stack_pacifico[] 的行将与 elevation_pacifico[] 的元素对应的单元格相同。
  • 感谢您的回答 jbaums。我还有一个疑问。很清楚如何使用raster::extract 和 SpatialPolygon 对象提取栅格数据,例如来自城市或国家的形状文件,但如何提取栅格每个像素的数据。我应该将栅格转换为 SpatialPolygon、Pixels 或 Points 对象,然后 extract 使用这些对象作为 y 的数据吗?
  • 如果你想要所有单元格的数据,那么下面的任何一个:as.data.frame(r); r[]; values(r)extract(r, seq_len(ncell(r))),其中rRaster* 对象,stack_pacificoelevation_pacifico。对于堆栈,这将为您提供data.frame,每层一列,每个单元格一行,对于单个栅格层,它将为您提供一个矢量。
  • 嗯,效果很好,我只将data.frame 选项添加到extract 函数中:extract(r, seq_len(ncell(r)), df=T) 给我一个data.frame,每个像素(或网格)都有一个ID 变量.我仍然想确定所有栅格中像素的识别。有没有办法获得像素的lan-long?我正在考虑尝试rasterToPoints 函数,但我也知道这个过程可能导致错误或丢失数据。有任何想法吗? (感谢您的广泛帮助!)
  • 所有单元格的坐标:coordinates(r)。对于特定单元格,xyFromCell(r, 1:10).

标签: r merge dataframe extract raster


【解决方案1】:

如果两个栅格对象的分辨率、范围和坐标系相同,则单元格将完美重叠。您可以通过查看坐标来确认这一点:

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 是单个栅格图层,它将有一个列;相比之下,替代方案将返回如果与单个栅格图层一起使用,则为简单矢量。)

as.data.frame(stack_pacifico) 的行与as.data.frame(elevation_pacifico) (or, equivalently, the elements ofvalues(elevation_pacifico)` 的行对应的单元格坐标相同。

【讨论】:

    【解决方案2】:

    或者这样做:

    s <- stack(elevation_pacifico, stack_pacifico)
    d <- values(s)
    

    【讨论】:

      猜你喜欢
      • 2022-06-15
      • 2021-09-09
      • 2020-07-19
      • 1970-01-01
      • 2021-10-17
      • 2018-09-15
      • 2020-04-06
      • 1970-01-01
      • 2018-09-17
      相关资源
      最近更新 更多