【问题标题】：R raster package split image into multiplesR光栅包将图像拆分为多个
【发布时间】：2015-06-29 08:37:16
【问题描述】：

我有一张如下图。它是 2579*2388 像素。让我们假设它的左下角是 0,0。从该图像中，我想创建多个图像，如下所示并将它们保存在工作文件夹中。每个图像的大小为 100*100 像素。每张图片将通过其左下角坐标保存。

第一张图片的左下角为 0,0。右上手角将在 100,100 并且图像将保存为 0-0.jpg
second 的左下角为 10,0。右上角角将在 110,100，图像将保存为 10-0.jpg
完成最后一行后，Y 坐标将移动 10。在第二行的情况下，第一张图像将在 0,10 并且该图像将保存为 0-10.jpg

最快的方法是什么？有没有可以非常快的 R 包？

我了解，对于当前图像，它会将其拆分为大约 257*238 个图像。但是我有足够的磁盘空间，我需要每张图片来执行文本检测。

【问题讨论】：

raster 包专用于“地理数据分析和建模”。
还有其他推荐的包吗？

标签： r image image-processing split raster

【解决方案1】：

这里使用“光栅”包的另一种方法。该函数对要切分的栅格进行空间聚合，聚合后的栅格像元转化为多边形，然后使用每个多边形的范围对输入栅格进行裁剪。

我确信有复杂而紧凑的方法可以做到这一点，但这种方法对我有用，而且我发现它也很直观。我希望你也觉得它有用。请注意，下面的第 4 部分和第 5 部分仅用于测试，它们不是函数的一部分。

第 1 部分：加载和绘制示例栅格数据

logo <- raster(system.file("external/rlogo.grd", package="raster"))
plot(logo,axes=F,legend=F,bty="n",box=FALSE)

第 2 部分：函数本身：

# The function spatially aggregates the original raster
# it turns each aggregated cell into a polygon
# then the extent of each polygon is used to crop
# the original raster.
# The function returns a list with all the pieces
# in case you want to keep them in the memory. 
# it saves and plots each piece
# The arguments are:
# raster = raster to be chopped            (raster object)
# ppside = pieces per side                 (integer)
# save   = write raster                    (TRUE or FALSE)
# plot   = do you want to plot the output? (TRUE or FALSE)
SplitRas <- function(raster,ppside,save,plot){
  h        <- ceiling(ncol(raster)/ppside)
  v        <- ceiling(nrow(raster)/ppside)
  agg      <- aggregate(raster,fact=c(h,v))
  agg[]    <- 1:ncell(agg)
  agg_poly <- rasterToPolygons(agg)
  names(agg_poly) <- "polis"
  r_list <- list()
  for(i in 1:ncell(agg)){
    e1          <- extent(agg_poly[agg_poly$polis==i,])
    r_list[[i]] <- crop(raster,e1)
  }
  if(save==T){
    for(i in 1:length(r_list)){
      writeRaster(r_list[[i]],filename=paste("SplitRas",i,sep=""),
                  format="GTiff",datatype="FLT4S",overwrite=TRUE)  
    }
  }
  if(plot==T){
    par(mfrow=c(ppside,ppside))
    for(i in 1:length(r_list)){
      plot(r_list[[i]],axes=F,legend=F,bty="n",box=FALSE)  
    }
  }
  return(r_list)
}

第 3 部分：测试功能

SplitRas(raster=logo,ppside=3,save=TRUE,plot=TRUE)
# in this example we chopped the raster in 3 pieces per side
# so 9 pieces in total
# now the raster pieces should be ready 
# to be processed in the default directory
# A feature I like about this function is that it plots
# the pieces in the original order.

第 4 部分：在每个部分上运行代码并将它们保存回目录中

# notice if you cropped a rasterbrick 
# use "brick" instead of "raster" to read
# the piece back in R
list2 <- list()
for(i in 1:9){ # change this 9 depending on your number of pieces
  rx <- raster(paste("SplitRas",i,".tif",sep=""))
  # piece_processed <- HERE YOU RUN YOUR CODE
  writeRaster(piece_processed,filename=paste("SplitRas",i,sep=""),
              format="GTiff",datatype="FLT4S",overwrite=TRUE)
}
# once a code has been ran on those pieces
# we save them back in the directory 
# with the same name for convenience

第 5 部分：让我们重新组合起来

# read each piece back in R
list2 <- list()
for(i in 1:9){ # change this 9 depending on your number of pieces
  rx <- raster(paste("SplitRas",i,".tif",sep=""))
  list2[[i]] <- rx
}
# mosaic them, plot mosaic & save output
list2$fun   <- max
rast.mosaic <- do.call(mosaic,list2)
plot(rast.mosaic,axes=F,legend=F,bty="n",box=FALSE)
writeRaster(rast.mosaic,filename=paste("Mosaicked_ras",sep=""),
            format="GTiff",datatype="FLT4S",overwrite=TRUE)

【讨论】：

一个非常优雅的解决方案。我想要同样的，但我正在使用具有 9 个 RasterLayers 的栅格堆栈。你会在代码中修改什么来保持堆叠的部分？
对不起Jecogeo，我不确定我是否理解这个问题。如果参数raster 填充了光栅堆栈，则输出片段也将保持为光栅堆栈。您是否使用 brick 读取了输入栅格堆栈？
是的，谢谢，@Shepherd。正如您所说，当作为光栅堆栈加载时，输出片段将保持为光栅堆栈。太好了！

【解决方案2】：

这是一种方法，通过 gdalUtils 使用 GDAL，如果需要，可以并行化。

library(gdalUtils)

# Get the dimensions of the jpg    
dims <- as.numeric(
  strsplit(gsub('Size is|\\s+', '', grep('Size is', gdalinfo('R1fqE.jpg'), value=TRUE)), 
           ',')[[1]]
)

# Set the window increment, width and height
incr <- 10
win_width <- 100
win_height <- 100

# Create a data.frame containing coordinates of the lower-left
#  corners of the windows, and the corresponding output filenames.
xy <- setNames(expand.grid(seq(0, dims[1], incr), seq(dims[2], 0, -incr)), 
               c('llx', 'lly'))
xy$nm <- paste0(xy$llx, '-', dims[2] - xy$lly, '.png')

# Create a function to split the raster using gdalUtils::gdal_translate
split_rast <- function(infile, outfile, llx, lly, win_width, win_height) {
  library(gdalUtils)
  gdal_translate(infile, outfile, 
                 srcwin=c(llx, lly - win_height, win_width, win_height))
}

将函数应用于单个窗口的示例：

split_rast('R1fqE.jpg', xy$nm[1], xy$llx[1], xy$lly[1], 100, 100)

将其应用于前 10 个窗口的示例：

mapply(split_rast, 'R1fqE.jpg', xy$nm[1:10], xy$llx[1:10], xy$lly[1:10], 100, 100)

使用 parLapply 并行运行的示例：

library(parallel)
cl <- makeCluster(4) # e.g. use 4 cores
clusterExport(cl, c('split_rast', 'xy')) 

system.time({
  parLapply(cl, seq_len(nrow(xy)), function(i) {
    split_rast('R1fqE.jpg', xy$nm[i], xy$llx[i], xy$lly[i], 100, 100)  
  })
})
stopCluster(cl)

【讨论】：

【解决方案3】：

这有点晚了，但可能对遇到此问题的其他人有用。 SpaDES 包有一个方便的函数，称为 splitRaster()，它可以满足您的需求。

一个例子：

library(raster)
library(SpaDES)

# Create grid
the_grid=raster(xmn=0, xmx=100, ymn=0, ymx=100, resolution=1)

# Set some values
the_grid[0:50,0:50] <- 1
the_grid[51:100,51:100] <- 2
the_grid[51:100,0:50] <- 3
the_grid[0:50,51:100] <- 4

这给了你这个：现在使用SpaDES package 进行拆分。根据您想要沿 x 和 y 轴的图块数量设置 nx 和 ny - 如果我们想要 4 个图块，请将它们设置为 nx=2 和 ny=2。如果您没有设置path，它应该将文件写入您的当前目录。还有其他的东西，比如缓冲——见?splitRaster：

# Split into sections - saves automatically to path
sections=splitRaster(the_grid, nx=2, ny=2, path="/your_output_path/")

变量sections 是一个栅格列表，the_grid 的每个部分都有一个 - 访问它们：

split_1=sections[[1]]

如果您想专门保存它们，只需使用writeRaster()。

要再次创建组合栅格，请使用mergeRaster()。

【讨论】：

这看起来很棒。不幸的是，我将 150MB 的栅格分成 30 个部分，最后我得到了 35GB 的数据...
@TWest 这不是 R 解决方案，但请看这里的问答 - 它应该更好地处理内存（虽然我还没有测试过）：gis.stackexchange.com/questions/14712/…

【解决方案4】：

你可以使用gdal和r，如link所示。

然后您将修改第 23 行以进行适当的偏移以允许生成的图块之间重叠。

【讨论】：

这似乎令人困惑。另外我不确定它是否告诉如何保存这些图像......是否可以提供简单的示例或工作代码？

【解决方案5】：

没有找到专门使用r 的直接实现我使用以下方法使用raster 操作，这可能对其他人来说很有趣。它生成范围并将原始栅格裁剪为它们。希望这会有所帮助！

## create dummy raster
n <- 50
r <- raster(ncol=n, nrow=n, xmn=4, xmx=10, ymn=52, ymx=54)
projection(r) <- "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"
values(r)     <- 1:n^2+rnorm(n^2)


n.side <-  2  # number of tiles per side
dx     <- (extent(r)[2]- extent(r)[1])/ n.side  # extent of one tile in x direction
dy     <- (extent(r)[4]- extent(r)[3])/ n.side  # extent of one tile in y direction
xs     <- seq(extent(r)[1], by= dx, length= n.side) #lower left x-coordinates
ys     <- seq(extent(r)[3], by= dy, length= n.side) #lower left y-coordinates
cS     <- expand.grid(x= xs, y= ys)

## loop over extents and crop
for(i in 1:nrow(cS)) {
  ex1 <- c(cS[i,1], cS[i,1]+dx, cS[i,2], cS[i,2]+dy)  # create extents for cropping raster
  cl1 <- crop(r, ex1) # crop raster by extent
  writeRaster(x = cl1, filename=paste("test",i,".tif", sep=""), format="GTiff", overwrite=T) # write to file
}

## check functionality...
test <- raster(paste("test1.tif", sep=""))
plot(test)

【讨论】：

【解决方案6】：

我认为您需要为您的处理部分创建一个函数（我们称之为“fnc”）和一个列出您制作的瓷砖数量的表格（我们称之为“tile.tbl”）以及假设您的地理大数据名为“obj”

    obj=GDALinfo("/pathtodata.tif")
    tile.tbl <- getSpatialTiles(obj, block.x= your size of interest, return.SpatialPolygons=FALSE)

然后使用降雪包将其并行化。这是一个例子：

    library(snowfall)
    sfInit(parallel=TRUE, cpus=parallel::detectCores())
    sfExport("tile.tbl", "fnc")
    sfLibrary(rgdal)
    sfLibrary(raster)
    out.lst <- sfClusterApplyLB(1:nrow(tile.tbl), function(x){ fnc(x, tile.tbl) })
    sfStop()

详细解释请见HERE

【讨论】：

【解决方案7】：

@Shepherd 提出的方法可行，但是您需要事先定义图像将被划分成的补丁数量（ppside 参数）。我们经常需要获得相同大小的补丁（例如训练 CNN）。此外，将@Shepherd 方法应用于具有不同尺寸的图像列表可能会出现问题。

这是一种使用raster 包将图像划分为相等块的方法。

我们将创建一些玩具数据。 35x64 栅格。我们希望将栅格划分为 10x10 块。 patch_size = 10.

my_img = raster(nrow=64, ncol=35)
values(my_img) = 1:ncell(my_img)

我们定义了patchifyR 函数，它将执行以下操作：

patchifyR <- function(img, patch_size){
  # load raster 
  if (!require("raster")) install.packages("raster")
  suppressPackageStartupMessages({library(raster)})
  # create image divisible by the patch_size
  message(paste0("Cropping original image. ", "Making it divisible by ", patch_size, "."))
  x_max <- patch_size*trunc(nrow(img)/patch_size)
  y_max <- patch_size*trunc(ncol(img)/patch_size)
  img <- crop(img, extent(img, 1, x_max, 1, y_max))
  # initializers
  lx = 1; ly = 1; p = 1
  ls.patches <- list()
  ls.coordinates <- list()
  # extract patches
  for(i in 1:(nrow(img)/patch_size)){
    for(j in 1:(ncol(img)/patch_size)){
      ls.patches[[p]] <- crop(img, extent(img, lx, (lx+patch_size)-1, ly, (ly+patch_size)-1))
      ls.coordinates[[p]] <- as.character(paste0("P_", p, "_X0_", lx, "_X1_", (lx+patch_size)-1, "_Y0_", ly, "_Y1_", (ly+patch_size)-1))
      message(paste0("Patch ", p, " created. Coordinates: ", "X0_", lx, "_X1_", (lx+patch_size)-1, "_Y0_", ly, "_Y1_", (ly+patch_size)-1))
      p = p + 1
      ly = ly + patch_size
    }
    ly = 1
    lx = lx + patch_size
  }
  # merge results: $patches and $names
  message("Matching results ... ")
  patchify <- list("patches"=ls.patches, "names"=ls.coordinates)
  # return
  message("Successfully completed.")
  return(patchify)
}

首先，我们将剪切光栅以创建可被patch_size 整除的图像。如果我们的图块的重叠大于patch_size，我们将不会丢失信息。

然后我们使用crop 函数和双for 循环来循环坐标来创建不同的补丁。

patchifyR 函数接收两个参数： img： 输入图像。 RasterLayer 或 RasterStack。 patch_size： 路径大小。 10 创建 10x10 像素的补丁。

patchifyR 函数返回两个列表： patchify$patches： 包含所有补丁的列表（RasterLayer 或 RasterStack） patchify$names： 包含每个补丁的名称及其坐标。

让我们测试一下这个功能：

my_patches <- patchifyR(img=my_img, patch_size=10)

最后我们可以使用以下命令将结果保存到磁盘：

output_directory <- "D:/my_output/"

for(i in 1:length(my_patches$patches)){
  writeRaster(my_patches$patches[[i]], paste0(output_directory, my_patches$names[[i]], ".tif"), drivername="Gtiff", overwrite=TRUE)
}

【讨论】：