Разделите RGB растр на несколько частей и сохраните файл в GeoTiff с условиями - PullRequest
Новая
       48

Разделите RGB растр на несколько частей и сохраните файл в GeoTiff с условиями

0 голосов
/ 08 июля 2019

Я хотел бы разбить RGB-растр на несколько частей и сохранить каждое изображение в формате GeoTiff, но в имени файла мне нужно количество точек (pts) внутри каждого фрагмента. Мой код работает отлично, когда я использовал одно растровое изображение, но со стеком или кирпичным объектом не работает. В качестве окончательного вывода в моей директории мне нужно девять файлов GeoTiff (я разделил исходный растр на 9 частей) с именами файлов: SplitRas1points0.tif ... SplitRas9points0.tif

Мой код: 

#Packages
library(raster)
library(sp)
library(rgdal)

## Create a simulated RBG rasters
r <- raster(nc=30, nr=30)
r <- setValues(r, round(runif(ncell(r))* 255))
g <- raster(nc=30, nr=30)
g <- setValues(r, round(runif(ncell(r))* 255))
b <- raster(nc=30, nr=30)
b <- setValues(r, round(runif(ncell(r))* 255))
rgb<-stack(r,g,b)
plotRGB(rgb,
        r = 1, g = 2, b = 3)

##Given interesting points coordinates
xd     <- c(-24.99270,45.12069,99.40321,73.64419)
yd  <- c(-45.435267,-88.369745,-7.086949,44.174530)
pts <- data.frame(xd,yd)
pts_s<- SpatialPoints(pts)
points(pts_s, col="black", pch=16)


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


#Slip RGB raster in 3 parts
splitRBG<-SplitRas(raster=rgb,ppside=3,save=FALSE,plot=FALSE)

Здесь ОК, есть 9 растров !! 

[[1]]
class      : RasterBrick 
dimensions : 10, 10, 100, 3  (nrow, ncol, ncell, nlayers)
resolution : 12, 6  (x, y)
extent     : -180, -60, 30, 90  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
source     : memory
names      : layer.1, layer.2, layer.3 
min values :       0,       0,       2 
max values :     249,     253,     252 

...

[[9]]
class      : RasterBrick 
dimensions : 10, 10, 100, 3  (nrow, ncol, ncell, nlayers)
resolution : 12, 6  (x, y)
extent     : 60, 180, -90, -30  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
source     : memory
names      : layer.1, layer.2, layer.3 
min values :       3,       3,       3 
max values :     254,     251,     248  



#Count number of points inside each piece of raster
res<-NULL
for(i in 1:3){
  pointcount = function(r, pts){
  # make a raster of zeroes like the input
  r2 = r
  r2[] = 0
  # get the cell index for each point and make a table:
  counts = table(cellFromXY(r,pts))
  # fill in the raster with the counts from the cell index:
  r2[as.numeric(names(counts))] = counts
  return(r2)
 }
r2 = pointcount(splitRBG[i], pts_s)
res<-rbind(res,c(r2))
}
#

Первая проблема, мой объект стека рассматривает список, и у меня есть вывод: 

Error in (function (classes, fdef, mtable)  : 
  unable to find an inherited method for function ‘cellFromXY’ for signature ‘"list"’


#Run a code on each piece with the number of points inside each piece in the file name (res[i])
list2 <- list()
for(i in 1:3){ # 3 is the number of pieces
  rx <- raster(paste("splitRBG",i,".tif",sep=""))
writeRaster(splitRBG,filename=paste("SplitRas",i,"points",res[i],sep=""),
              format="GTiff",datatype="FLT4S",overwrite=TRUE)
}
#

Вторая проблема. Хорошо, если у меня нет точечных чисел, являющихся причиной вышеуказанной ошибки, но мой объект splitRBG не распознается как RasterLayer, посмотрите мой вывод: 

Error in .rasterObjectFromFile(x, band = band, objecttype = "RasterLayer",  : 
  Cannot create a RasterLayer object from this file. (file does not exist)

Есть идеи, пожалуйста?

1 Ответ

0 голосов
/ 26 июля 2019

Теперь у меня есть rcount<-length(intersect(SpatialPoints(pts_s), r_list[[i]])) внутри первой функции, и проблема решена !! 

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)){
      rcount<-length(intersect(SpatialPoints(pts_s), r_list[[i]])) #Count number of points inside each piece of raster
      writeRaster(r_list[[i]],filename=paste(rcount,"_sample_",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)
}


#Slip RGB raster in 3 parts
splitRBG<-SplitRas(raster=rgb,ppside=3,save=TRUE,plot=FALSE)
#
...