1

I've like to split a RGB raster in multiple pieces and save each image in GeoTiff format, but in the file name I need the count of points (pts) inside each piece. My code runs perfectly when I used a single raster image, but with stack or brick object doesn't work. I've like as final output in my directory nine GeoTiff files (I divided the original raster in 9 parts) with files names: SplitRas1points0.tif ... SplitRas9points0.tif

My code is:

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

Here OK I have 9 raster!!

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

First problem, my stack object is considering a list and I have the output:

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

Second problem, Ok if I don't have the numbers of point a cause of the error above, but my splitRBG object is don't recognise as RasterLayer, look my output:

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

Any ideas please?

  • Do you get an error message? No output? Computer catch fire? Why are you creating rx and then doing nothing with it? Do you mean to save this raster? – Spacedman Jul 8 at 15:20
  • Spacedman, I make all the comments about errors in the question. Thanks! – Leprechault Jul 9 at 13:30
0

Now, I have rcount<-length(intersect(SpatialPoints(pts_s), r_list[[i]])) inside the first function and the problem is solved!!

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.