如何更改 R 中栅格图层的分辨率
How to change the resolution of a raster layer in R
我正在使用 R 中的多个高分辨率栅格图层。对于我 运行 的某些分析来说，详细程度过高，因此我想通过降低分辨率来加快速度。
坐标系是UTM所以单位是米。分辨率说是30, 30(x, y)。所以这里的分辨率好像是30m。
有人可以告诉我如何将分辨率改为 120m 吗？我已经阅读了 resample() 和 projectRaster() 函数的帮助，但它们似乎需要具有所需分辨率的模板栅格，而我没有。
这是我的栅格图层之一的示例：
alt.utm
class : RasterLayer
dimensions : 4572, 2495, 11407140 (nrow, ncol, ncell)
resolution : 30, 30 (x, y)
extent : 421661, 496511, 4402939, 4540099 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=13 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
data source : in memory
names : layer
values : 1485.127, 4275.202 (min, max)
这是一个如何执行此操作的示例。 (link to original)
#########################################################################
# SubsampleImageRaster.r
#
# This function demonstrates resampling of raster images to a new
# spatial resolution using the R raster package.
#
# Author: Rick Reeves
# Date created: 6-October- 2010
# Date modified:
# NCEAS
#
#########################################################################
#
SubsampleImageRaster <- function()
{
library(raster)
sOutFile <- ""
resampleFactor <- 4 # For test, subsample incoming image by factor of 10
# Read the mosaic components, stored in a subfolder, into a raster object list.
# Within the same loop, obtain the (geographic) extent of each component.
# Note: these images do not have same spatial extent, so they cant be stored
# in a rasterStack. Instead, use a list of rasterLayers.
setwd("./ForUseCase")
inFiles <- list.files(pattern="*.tif")
nFiles <- length(inFiles)
inputRaster <- raster()
CoarseResampRaster <- raster()
FineResampRaster <- raster()
for (iCtr in 1 : nFiles)
{
message(sprintf("resampling file: %s",inFiles[iCtr]))
inputRaster <- raster(inFiles[iCtr])
# The aggregate() / disaggregate methods resample rasters to COARSER (bigger cells)
# and FINER (smaller cells) resolutions, respectively
CoarseResampRaster <- aggregate(inputRaster,fact=resampleFactor,fun=mean)
sOutFile <- sprintf("CoarseSubsamp%s",inFiles[iCtr])
writeRaster(CoarseResampRaster,filename=sOutFile,format="GTiff",datatype="INT1U",overwrite=TRUE)
FineResampRaster <- disaggregate(inputRaster,fact=resampleFactor,fun=mean)
sOutFile <- sprintf("FineSubsamp%s",inFiles[iCtr])
writeRaster(FineResampRaster,filename=sOutFile,format="GTiff",datatype="INT1U",overwrite=TRUE)
}
message("resample demo")
browser()
# second method: use the resample() method from raster package
# Simple example:
# This code produces a resampled raster, 's',
# with correct resampled values. e.g.;
# s[] prints a vector of resampled cell values.
r <- raster(nrow=3, ncol=3)
r[] <- 1:ncell(r)
s <- raster(nrow=10, ncol=10)
s <- resample(r, s, method='bilinear')
# Useful example:
# Resample a satellite image, stored in a GeoTiff file
# into a NEW raster with 2x spatial resolution in
# both dimensions (four times the number of cells)
# Here is the technique:
# 1) Create a new raster object with the correct 'resampled' number of cells.
# 2) Set the extent (geographic 'bounding box') of the new raster
# to the extent of the original raster
# 3) Generate the resampled raster.
resampleFactor <- .5 # reduce the cell size by 50% and double the number of rows and columns.
inputRaster <- raster("TmB50MosaicImg1.tif")
inCols <- ncol(inputRaster)
inRows <- nrow(inputRaster)
resampledRaster <- raster(ncol=(inCols / resampleFactor), nrow=(inRows / resampleFactor))
extent(resampledRaster) <- extent(inputRaster)
# The resample method will write the resampled raster image to a NEW disk file..
resampledRaster <- resample(inputRaster,resampledRaster,datatype="INT1U",method='bilinear',filename="testOutResamp.tif",overwrite=TRUE)
# Or, use writeRaster method to create the output file.
writeRaster(resampledRaster,filename="ResampleProduct.tif",format="GTiff",datatype="INT1U",overwrite=TRUE)
# end
}
您可以使用 aggregate or disaggregate.
library(raster)
#get some sample data
data(meuse.grid)
gridded(meuse.grid) <- ~x+y
meuse.raster <- raster(meuse.grid)
res(meuse.raster)
#[1] 40 40
#aggregate from 40x40 resolution to 120x120 (factor = 3)
meuse.raster.aggregate <- aggregate(meuse.raster, fact=3)
res(meuse.raster.aggregate)
#[1] 120 120
#disaggregate from 40x40 resolution to 10x10 (factor = 4)
meuse.raster.disaggregate <- disaggregate(meuse.raster, fact=4)
res(meuse.raster.disaggregate)
#[1] 10 10
我最近需要降低 ggmap 对象的分辨率。这涉及提取和转换 ggmap 栅格（使用 Robin Lovelace 的 ggmap_rast()），按照本线程中讨论的那样聚合栅格，然后用下面的较低分辨率栅格替换 ggmap 高分辨率栅格。希望这有用：
original_map <- get_map("New York", scale = 1) #download a map at lowest resolution
#ggmap_rast function courtesy of Robin Lovelace
original_map.rast1 <- ggmap_rast(original_map) #extract RasterStack from ggmap object
original_map.rast2 <- aggregate(original_map.rast, 2) #compress raster
rast2_length <- sqrt(length(original_map.rast2)/3) ## find the number of cells in compressed raster
map <- ggmap::ggmap(original_map)
# from
r <- map$layers[[2]]$geom_params$raster #pull hex raster pixel values from ggmap object
#x <- r[,] #assign values to a variable (probably unnecessary)
xx <- r[1:rast2_length,1:rast2_length] #filler values for raster to be created
rgb_map <- original_map.rast2
for(i in 1:rast2_length){
for(j in 1:rast2_length){
k=i*rast2_length+j
#many rgv values are non-integers; rgb2hex requires integer
red = as.integer(round(rgb_map$layer.1[k],0))
green = as.integer(round(rgb_map$layer.2[k],0))
blue = as.integer(round(rgb_map$layer.3[k],0))
#rgb2hex from ggtern package
xx[i,j] <- rgb2hex(red,green,blue)
#rgb2hex is slow; export, calculate in excel, import is faster, sadly
#xx[i,j] <- as.character(rgb_from_excel$V1[k])
}
}
map$layers[[2]]$geom_params$raster <- xx
map
我已经尝试了 3 种不同的选项来放大 DEM 文件。首先我像这样使用 gdal_translate：
from CMD:
gdal_translate -tr 0.1 0.1 "C:\dem.tif" "C:\dem_0.1.tif"
在 R 中：
res(raster('C:\dem_0.1.tif')
[1] 0.1 0.1
然后我在 R:
中尝试了 raster 包中的 aggregate 和 resample 命令
#
r <- raster("dem.tif")
> res(r)
[1] 0.0002777778 0.0002777778 # original dem resolution
#
r_up_0.1 <- aggregate(r, fact = 0.1/res(r)) # aggregate output
> res(r_up_0.1)
[1] 0.1 0.1
#
s <- raster(nrow = 10, ncol = 10)
extent(s) <- extent(r)
s <- resample(r, s, method = 'bilinear') # resample output
> res(s)
[1] 0.1000278 0.1000278
这些是输出：
3 个输出被放大到 res = 0.1x0.1，它们之间存在一些差异，但我打算使用 gdal 输出。
希望对您有所帮助。
这是一个迟到的答案，但创建模板栅格并使用 projectRaster 非常简单。您可以将模板光栅分辨率设置为您想要的任何值，如果您想要矩形而不是方形单元格，可以包括两个值：
# Create the template raster, with 120m cells
TEMPLATE.RASTER <- raster(extent(ORIGINAL.RASTER), resolution = 120,
crs = st_crs(ORIGINAL.RASTER)$proj4string)
# Project the original raster to the new resolution
PROJECTED.RASTER <- projectRaster(from = ORIGINAL.RASTER, to = TEMPLATE.RASTER)
我喜欢创建一个模板，因为它也会复制目标 CRS。但是，您实际上可以跳过创建模板光栅并直接在 projectRaster 中设置分辨率和 CRS。请务必查看帮助文件，因为有计算方法等选项可能很重要。
这些天我们可以使用 terra，替代 raster 包
示例数据
library(terra)
f <- system.file("ex/elev.tif", package="terra")
r <- rast(f)
聚合栅格像元。从 30 米到 120 米是 4 倍
a1 <- aggregate(r, 4)
您可以对行和列使用不同的因子，也可以使用不同的聚合函数（默认为“均值”）
a2 <- aggregate(r, c(2,3), fun=sum)
你也可以反其道而行之：
d <- disagg(a2, 2)
聚合只能合并整个单元格。但是要合并来自不同来源的栅格数据，您可能需要匹配未对齐的栅格的几何形状。在这种情况下，您可以使用 resample
非对齐示例 SpatRaster:
x <- rast(r)
res(x) <- 0.01
解决方案
x <- resample(r, x)
也可能是您想将栅格数据转换为具有另一个坐标参考系统（“地图投影”）的几何图形。你可以这样做：
u1 <- project(r, "+proj=utm +zone=32")
但是，与矢量数据不同，栅格数据的转换定义不明确。因此，首选方法是提供一个您希望输出与之对齐的模板。这通常是您已经从另一个数据源获得的 SpatRaster。但在这里我为示例创建了一个：
temp <- rast(xmin=264000, xmax=324000, ymin=5479000, ymax=5565000, res=100, crs="+proj=utm +zone=32")
现在使用它
u2 <- project(r, temp)
进一步reading
我正在使用 R 中的多个高分辨率栅格图层。对于我 运行 的某些分析来说，详细程度过高，因此我想通过降低分辨率来加快速度。
坐标系是UTM所以单位是米。分辨率说是30, 30(x, y)。所以这里的分辨率好像是30m。
有人可以告诉我如何将分辨率改为 120m 吗？我已经阅读了 resample() 和 projectRaster() 函数的帮助，但它们似乎需要具有所需分辨率的模板栅格，而我没有。
这是我的栅格图层之一的示例：
alt.utm
class : RasterLayer
dimensions : 4572, 2495, 11407140 (nrow, ncol, ncell)
resolution : 30, 30 (x, y)
extent : 421661, 496511, 4402939, 4540099 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=13 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
data source : in memory
names : layer
values : 1485.127, 4275.202 (min, max)
这是一个如何执行此操作的示例。 (link to original)
#########################################################################
# SubsampleImageRaster.r
#
# This function demonstrates resampling of raster images to a new
# spatial resolution using the R raster package.
#
# Author: Rick Reeves
# Date created: 6-October- 2010
# Date modified:
# NCEAS
#
#########################################################################
#
SubsampleImageRaster <- function()
{
library(raster)
sOutFile <- ""
resampleFactor <- 4 # For test, subsample incoming image by factor of 10
# Read the mosaic components, stored in a subfolder, into a raster object list.
# Within the same loop, obtain the (geographic) extent of each component.
# Note: these images do not have same spatial extent, so they cant be stored
# in a rasterStack. Instead, use a list of rasterLayers.
setwd("./ForUseCase")
inFiles <- list.files(pattern="*.tif")
nFiles <- length(inFiles)
inputRaster <- raster()
CoarseResampRaster <- raster()
FineResampRaster <- raster()
for (iCtr in 1 : nFiles)
{
message(sprintf("resampling file: %s",inFiles[iCtr]))
inputRaster <- raster(inFiles[iCtr])
# The aggregate() / disaggregate methods resample rasters to COARSER (bigger cells)
# and FINER (smaller cells) resolutions, respectively
CoarseResampRaster <- aggregate(inputRaster,fact=resampleFactor,fun=mean)
sOutFile <- sprintf("CoarseSubsamp%s",inFiles[iCtr])
writeRaster(CoarseResampRaster,filename=sOutFile,format="GTiff",datatype="INT1U",overwrite=TRUE)
FineResampRaster <- disaggregate(inputRaster,fact=resampleFactor,fun=mean)
sOutFile <- sprintf("FineSubsamp%s",inFiles[iCtr])
writeRaster(FineResampRaster,filename=sOutFile,format="GTiff",datatype="INT1U",overwrite=TRUE)
}
message("resample demo")
browser()
# second method: use the resample() method from raster package
# Simple example:
# This code produces a resampled raster, 's',
# with correct resampled values. e.g.;
# s[] prints a vector of resampled cell values.
r <- raster(nrow=3, ncol=3)
r[] <- 1:ncell(r)
s <- raster(nrow=10, ncol=10)
s <- resample(r, s, method='bilinear')
# Useful example:
# Resample a satellite image, stored in a GeoTiff file
# into a NEW raster with 2x spatial resolution in
# both dimensions (four times the number of cells)
# Here is the technique:
# 1) Create a new raster object with the correct 'resampled' number of cells.
# 2) Set the extent (geographic 'bounding box') of the new raster
# to the extent of the original raster
# 3) Generate the resampled raster.
resampleFactor <- .5 # reduce the cell size by 50% and double the number of rows and columns.
inputRaster <- raster("TmB50MosaicImg1.tif")
inCols <- ncol(inputRaster)
inRows <- nrow(inputRaster)
resampledRaster <- raster(ncol=(inCols / resampleFactor), nrow=(inRows / resampleFactor))
extent(resampledRaster) <- extent(inputRaster)
# The resample method will write the resampled raster image to a NEW disk file..
resampledRaster <- resample(inputRaster,resampledRaster,datatype="INT1U",method='bilinear',filename="testOutResamp.tif",overwrite=TRUE)
# Or, use writeRaster method to create the output file.
writeRaster(resampledRaster,filename="ResampleProduct.tif",format="GTiff",datatype="INT1U",overwrite=TRUE)
# end
}
您可以使用 aggregate or disaggregate.
library(raster)
#get some sample data
data(meuse.grid)
gridded(meuse.grid) <- ~x+y
meuse.raster <- raster(meuse.grid)
res(meuse.raster)
#[1] 40 40
#aggregate from 40x40 resolution to 120x120 (factor = 3)
meuse.raster.aggregate <- aggregate(meuse.raster, fact=3)
res(meuse.raster.aggregate)
#[1] 120 120
#disaggregate from 40x40 resolution to 10x10 (factor = 4)
meuse.raster.disaggregate <- disaggregate(meuse.raster, fact=4)
res(meuse.raster.disaggregate)
#[1] 10 10
我最近需要降低 ggmap 对象的分辨率。这涉及提取和转换 ggmap 栅格（使用 Robin Lovelace 的 ggmap_rast()），按照本线程中讨论的那样聚合栅格，然后用下面的较低分辨率栅格替换 ggmap 高分辨率栅格。希望这有用：
original_map <- get_map("New York", scale = 1) #download a map at lowest resolution
#ggmap_rast function courtesy of Robin Lovelace
original_map.rast1 <- ggmap_rast(original_map) #extract RasterStack from ggmap object
original_map.rast2 <- aggregate(original_map.rast, 2) #compress raster
rast2_length <- sqrt(length(original_map.rast2)/3) ## find the number of cells in compressed raster
map <- ggmap::ggmap(original_map)
# from
r <- map$layers[[2]]$geom_params$raster #pull hex raster pixel values from ggmap object
#x <- r[,] #assign values to a variable (probably unnecessary)
xx <- r[1:rast2_length,1:rast2_length] #filler values for raster to be created
rgb_map <- original_map.rast2
for(i in 1:rast2_length){
for(j in 1:rast2_length){
k=i*rast2_length+j
#many rgv values are non-integers; rgb2hex requires integer
red = as.integer(round(rgb_map$layer.1[k],0))
green = as.integer(round(rgb_map$layer.2[k],0))
blue = as.integer(round(rgb_map$layer.3[k],0))
#rgb2hex from ggtern package
xx[i,j] <- rgb2hex(red,green,blue)
#rgb2hex is slow; export, calculate in excel, import is faster, sadly
#xx[i,j] <- as.character(rgb_from_excel$V1[k])
}
}
map$layers[[2]]$geom_params$raster <- xx
map
我已经尝试了 3 种不同的选项来放大 DEM 文件。首先我像这样使用 gdal_translate：
from CMD:
gdal_translate -tr 0.1 0.1 "C:\dem.tif" "C:\dem_0.1.tif"
在 R 中：
res(raster('C:\dem_0.1.tif')
[1] 0.1 0.1
然后我在 R:
中尝试了raster 包中的 aggregate 和 resample 命令
#
r <- raster("dem.tif")
> res(r)
[1] 0.0002777778 0.0002777778 # original dem resolution
#
r_up_0.1 <- aggregate(r, fact = 0.1/res(r)) # aggregate output
> res(r_up_0.1)
[1] 0.1 0.1
#
s <- raster(nrow = 10, ncol = 10)
extent(s) <- extent(r)
s <- resample(r, s, method = 'bilinear') # resample output
> res(s)
[1] 0.1000278 0.1000278
这些是输出：
res = 0.1x0.1，它们之间存在一些差异，但我打算使用 gdal 输出。
希望对您有所帮助。
这是一个迟到的答案，但创建模板栅格并使用 projectRaster 非常简单。您可以将模板光栅分辨率设置为您想要的任何值，如果您想要矩形而不是方形单元格，可以包括两个值：
# Create the template raster, with 120m cells
TEMPLATE.RASTER <- raster(extent(ORIGINAL.RASTER), resolution = 120,
crs = st_crs(ORIGINAL.RASTER)$proj4string)
# Project the original raster to the new resolution
PROJECTED.RASTER <- projectRaster(from = ORIGINAL.RASTER, to = TEMPLATE.RASTER)
我喜欢创建一个模板，因为它也会复制目标 CRS。但是，您实际上可以跳过创建模板光栅并直接在 projectRaster 中设置分辨率和 CRS。请务必查看帮助文件，因为有计算方法等选项可能很重要。
这些天我们可以使用 terra，替代 raster 包
示例数据
library(terra)
f <- system.file("ex/elev.tif", package="terra")
r <- rast(f)
聚合栅格像元。从 30 米到 120 米是 4 倍
a1 <- aggregate(r, 4)
您可以对行和列使用不同的因子，也可以使用不同的聚合函数（默认为“均值”）
a2 <- aggregate(r, c(2,3), fun=sum)
你也可以反其道而行之：
d <- disagg(a2, 2)
聚合只能合并整个单元格。但是要合并来自不同来源的栅格数据，您可能需要匹配未对齐的栅格的几何形状。在这种情况下，您可以使用 resample
非对齐示例 SpatRaster:
x <- rast(r)
res(x) <- 0.01
解决方案
x <- resample(r, x)
也可能是您想将栅格数据转换为具有另一个坐标参考系统（“地图投影”）的几何图形。你可以这样做：
u1 <- project(r, "+proj=utm +zone=32")
但是，与矢量数据不同，栅格数据的转换定义不明确。因此，首选方法是提供一个您希望输出与之对齐的模板。这通常是您已经从另一个数据源获得的 SpatRaster。但在这里我为示例创建了一个：
temp <- rast(xmin=264000, xmax=324000, ymin=5479000, ymax=5565000, res=100, crs="+proj=utm +zone=32")
现在使用它
u2 <- project(r, temp)
进一步reading