26

我正在使用的 R 中有几个高分辨率栅格图层。对于我正在运行的某些分析，详细程度过高，因此我想通过降低分辨率来加快速度。

坐标系是 UTM，所以单位是米。分辨率说是 30, 30 (x, y)。所以看起来这里的分辨率是30m。

有人可以告诉我如何将分辨率更改为 120m 吗？我已阅读 resample() 和 projectRaster() 函数的帮助，但它们似乎需要具有所需分辨率的模板栅格，而我没有。

这是我的栅格图层之一的示例：

alt.utm
类：RasterLayer
尺寸：4572、2495、11407140（nrow、ncol、ncell）
分辨率：30、30（x、y）
范围：421661、496511、4402939、4540099（xmin、xmax、ymin、ymax）
坐标。参考。: +proj=utm +zone=13 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
数据源 : 在内存
名称中 : 层
值 : 1485.127, 4275.202 (min ， 最大限度）

4

6 回答 6

48

您可以使用聚合或分解。

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
于 2015-08-29T09:11:01.130 回答
10

我尝试了 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 中包中的aggregate和resample命令：raster

#
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 个输出被放大到 ares = 0.1x0.1并且它们之间存在一些差异，但我打算使用 gdal 输出。

希望这可以帮助。

于 2018-06-10T01:01:37.440 回答
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
}
于 2015-08-28T20:03:02.327 回答
4

这些天我们可以使用terra，更换raster包装

示例数据

library(terra)
f <- system.file("ex/elev.tif", package="terra")
r <- rast(f)

聚合栅格像元。从 30 到 120m 是 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)

进一步阅读

于 2021-11-25T18:27:12.553 回答
2

这是一个较晚的答案，但创建模板栅格并使用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. 请检查帮助文件，因为计算方法等选项可能很重要。

于 2021-11-25T10:17:24.180 回答
-1

我最近需要降低 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 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
于 2018-03-08T23:28:11.533 回答