如何在R中更改栅格图层的分辨率

Question

如何在R中更改栅格图层的分辨率

rresolutionrasterutm

36

我在R中有几个高分辨率的栅格图层，正在处理中。对于一些分析而言，细节级别太高了，我希望通过降低分辨率来加快速度。

坐标系统是UTM，单位为米。分辨率是30,30（x,y）。因此，分辨率似乎是30米。

请问有人可以建议我如何将分辨率改为120米？我已经阅读了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)

- NRP

6个回答

12

现在我们可以使用 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)

更多阅读

- Robert Hijmans

11

我已经尝试了三种不同的选项来升级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输出结果。

希望这有所帮助。

- noriega

三种方法的前后总数相比如何？我经常这样做，单元格值（和栅格总数）可能会有很大的变化。 - Sam

5

这里有一个做法示例。 (原文链接)

  #########################################################################
# 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
}

- scribbles

谢谢你的笔记，这对于在任何方向上操纵分辨率都很有用。 - NRP

3

这是一个晚回答，但创建模板栅格并使用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)

我喜欢创建一个模板，因为它也会复制目标坐标系。但是，您实际上可以跳过创建模板栅格并直接在 projectRaster 中设置分辨率和坐标系。请查看帮助文件，因为有一些计算方法等选项可能很重要。

- EcologyTom

-1

最近我有一个要求，需要降低ggmap对象的分辨率。这涉及到提取和转换ggmap栅格（使用Robin Lovelace的ggmap_rast()函数），按照这个讨论中所述进行栅格聚合，然后用下面的低分辨率栅格替换高分辨率栅格。希望这对您有用：

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 https://stackoverflow.com/questions/44225063/plot-ggmap-image-over-raster
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

- PFE

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- jhhwilliams · Accepted Answer

你可以使用 aggregate 或 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