【问题标题】:projectRaster results in different output when projecting with its own crs使用自己的 crs 进行投影时,projectRaster 会产生不同的输出
【发布时间】:2020-04-21 21:02:03
【问题描述】:

对于一个函数,我想使用用户设置的 CRS 重新投影光栅输入 - 使用裁剪和蒙版后的输出本身。我认为使用现有的 crs 进行投影不会做任何事情,只是返回输入栅格。令我惊讶的是,事实并非如此。下面是一个可重现的例子

创建虚拟栅格:

library(raster)

# Download country polygin, in this case Malawi
mwi <- raster::getData("GADM", country = "MWI", level = 0)

# Create dummy raster
grid <- raster::raster() # 1 degree raster
grid <- raster::disaggregate(grid, fact = 12) # 5 arcmin
grid <- raster::crop(grid, mwi)
values(grid) <- rep(1, ncell(grid)) # Set a value

# The input raster with dimensions 94,39,3666
grid <- raster::mask(grid, mwi)
plot(grid)
grid

# Reproject the raster using its own crs. I use ngb as it is a categorical variable.
# This raster has dimensions 102, 47, 4794 so it seems a lot of white space (NA) is added.
own_crs <- crs(grid)
grid_reproj <- raster::projectRaster(grid, crs = own_crs, method = "ngb")
plot(grid_reproj)
grid_reproj

# To remove the white space I use trim
# This results in a waster with dimensions 93, 39, 3627
grid_trim <- raster::trim(grid_reproj)
plot(grid_trim)
grid_trim

# I also decided to compare the maps visually with mapview
library(mapview)

# There seems to be a trim function in mapview which I set to FALSE
# Also use the browser for easy viewing
options(viewer = NULL)
mapviewOptions(trim = FALSE)
mapviewGetOption("trim")
mapview(grid, col.regions = "green", legend = F) +
  mapview(grid_reproj, col.regions = "red", legend = F) +
  mapview(grid_trim, col.regions = "blue", legend = F) 

比较地图我观察到两件事:

(1) grid 和 grid_trim 除了单个网格单元外几乎相同 在上面。也许这是由于四舍五入,

(2) grid_reproj 具有更大的维度和不同的范围。它也是 与其他地图相比,该地图似乎略有偏移。这 由修剪校正,所以我假设这些实际上是 NA 细胞,差异可能是 与 mapview 如何显示地图有关。

因此,我的主要问题是,当 rasterProject 使用 一样的程度?为什么结果会不同,即使在修剪之后?

【问题讨论】:

    标签: r gis raster map-projections


    【解决方案1】:

    要以正确的方式投影栅格数据,您必须提供另一个 Raster* 对象,而不仅仅是 crs。

    如果您提供 crs,则会计算一个新范围,该范围比预期的要大一些,以避免丢失数据。当然,在 crs 实际上相同的这种奇怪情况下,将输入保持不变返回是有意义的。

    【讨论】:

    • 感谢罗伯特的解释。这当然可以澄清。也许您可以更新 projectRaster 的文档以反映这一点,因为它当前显示:“您可以通过提供新投影作为单个参数来做到这一点,在这种情况下,函数设置新对象的范围和分辨率”。因此它建议人们只能使用 CRS 参数,事实上,我在教程中阅读过建议以这种方式使用它来简单地重新投影栅格。
    猜你喜欢
    • 1970-01-01
    • 2021-02-12
    • 2019-08-09
    • 2021-05-15
    • 2013-03-16
    • 2020-09-04
    • 2020-06-11
    • 2015-12-14
    • 1970-01-01
    相关资源
    最近更新 更多