1

I have a shapefile which contains arround 800 circular polygons with diameter of 60 meters each. Raster is 20 meter resolution so I would have multiple cell values extracted for each polygon. I am using an extract function on multiple raster layers in a for loop, and am putting the results in the dataframe. I know i can get the cellnumbers within the extract function, but I would also like to obtain the feature values (in my case a column in attribute table containing specific polygon code) directly from the shapefile for each pixel value extracted. Thus I would have in one row a pixel cell number, its extracted value and the polygon code.

my code:

dirs=list.dirs(full.names=TRUE)
dirs=dirs[grepl("R20",dirs)]
dirs=as.list(dirs)

DataToWrite=list()
sumirana=data.frame()
for (j in 1:length(dirs)) {
  setwd(paste("E:/",as.character(dirs[j]),sep="/"))
  files=list.files(full.names=FALSE,pattern="*20m*.tif$")
  for (i in 1:length(files)) {
      curRaster=raster(paste(getwd(),files[i],sep="/"))
      rasterOut=(na.omit(extract(curRaster,azo_plohe, cellnumbers=TRUE, 
      weights=TRUE, df=TRUE)))
      if((length(rasterOut) > 1)) {

         DataToWrite=c(DataToWrite,rasterOut[3])

         } else {

         lista=list(rep("x",15))
         tablica_folder=as.data.frame(t(unlist(lista)))      
         colnames(tablica_folder)=c(1:15)
         }
   }
   if((length(rasterOut) > 1)) {
      mylist=as.data.frame(DataToWrite)
      naziv_fajla=as.vector(rep((substr((names(DataToWrite)
      [1]),50,55)),nrow(mylist)))  
      tablica_folder=as.data.frame(c(rasterOut[2],naziv_fajla[1], 
      rasterOut[4],mylist),col.names=c(1:15))
      sumirana=rbind(sumirana,tablica_folder)
      DataToWrite=list()

      } else {
      sumirana=rbind(sumirana,tablica_folder)
      }
}
  • I think it may be more helpful (and you'll get better answers) if you provide the code for where you're stuck with your analysis & expected output. From your comments, do you need help with iterating through the folders? Or getting the raster values into a data.frame? – GISKid Nov 7 '17 at 21:34
  • HI! Thanx for the input. I do not need help with iterating or getting raster values to dataframe. This is a code that works well for me, I just want to upgrade it. I want to get the shapefile values (not raster) into the dataframe that is created at the end of my code (I think this is the place where I should put it : tablica_folder=as.data.frame(c(rasterOut[2],naziv_fajla[1], rasterOut[4],mylist),col.names=c(1:15))". However I do not know how to extract those shapefile values within my code. I want each extracted value from the shapefile polygon to correspond to the pixel extracted. – davor korman Nov 8 '17 at 11:16
  • For example, if the polygon that was used as a base for the pixel extraction has a code "CODE1" under the column "Polygon codes" which is a part of the shapefiles atribute table, I want all the pixels that were extracted (lets say the polygon has 30 meters diameter and has 7 pixels extracted) have the CODE1 added in my dataframe at the end of the code (which is sumirana = rbind(sumirana, tablica folder)). – davor korman Nov 8 '17 at 11:20
  • See my edits, you just need to merge your polygon to the new data.frame based on id. – GISKid Nov 9 '17 at 16:11
1

You have to unlist() the output of extract, and maintain that list-level grouping to know which object the value is from - this is a pain, and can be tricky for non-matches, so I put the workflow into package 'tabularaster' https://CRAN.R-project.org/package=tabularaster

You might find this is enough, where r is a raster and poly is a polygondataframe:

library(tabularaster)
cells <- cellnumbers(r, poly)
cells$value <- extract(r, cells$cell_)

So value is the pixel value, and cell_ is the pixel index, and object_ is the row-number of poly.

Then use object_ to get a value from the query layer:

## replace "polyID" with whatever column you want
cells$polyID <- poly$polyID[cells$object_]

For a case like this, making the cell_ explicit is overkill, but it works well for general extraction (i.e. time series coming in the door).

  • Thank you for your repply! I have updated the question with the code. I am not sure if I am able to implement this solution you proposed using the unlist option somewhere, because my whole code depends on keeping it in a list and then transforming it into a dataframe. – davor korman Nov 7 '17 at 13:21
1

This would be easier to answer with a reproducible set of data. I've created some here:

library(rgdal)
library(raster)
library(sp)
library(rgeos)
#Create fake raster data
r1 <- r2 <- r3 <- raster(res=20)
values(r1) <- runif(ncell(r1))
values(r2) <- runif(ncell(r2))
values(r3) <- runif(ncell(r3))
s <- stack(r1, r2, r3)
#name the rasters
names(s) <- c("raster1","vegetation", "raster3")
plot(s)
#set random sample for 10 points to be converted to polygons
set.seed(1)
rn<-runif(10, min=-150, max=150)
set.seed(999)
rn2<-runif(10, min=-150, max=150)
df<-data.frame(x=rn, y=rn2, id=seq(10))
cor<-df[,c(1,2)]
#convert to points
p<-SpatialPointsDataFrame(data=df, coords=cor)
#buffer points to create polygons
pb<-gBuffer(p, width=5, byid=TRUE)

plot(s[[1]])
plot(pb,add=T)

Once you have your rasters identified from your directory structure
the code below will take values from your polygons and put them into a data.frame with the name & value of the raster and the cell #. Since your rasters may be in different extents, stacking may not work but see this question "Handing multiple extent problem to create raster stack in r"

  #extract values from all rasters for each polygon
  #add new field with "code values"
  pb$code<-sample(10:20,10)
  ext_poly<-extract(s,pb,cellnumbers=TRUE, df=TRUE)
  ext_poly
#join polygon data to extracted values
  m<-merge(ext_poly,pb@data)

Here's the output with raster cell, raster cell values, and the joined polygon data (X,Y and the new field created above with "codes")

ID cell    raster1 vegetation   raster3          x         y code
1   1  114 0.47134842 0.11526912 0.7113075  -70.34740 -33.27859   16
2   1   96 0.71419069 0.64174118 0.8939626  -70.34740 -33.27859   16
3   2   62 0.56123275 0.97922882 0.6705206  -38.36283  24.91822   15
4   2   61 0.32030622 0.16661933 0.8966194  -38.36283  24.91822   15
5   5    5 0.48085049 0.84698577 0.5859292  -89.49542  86.02403   19
6   5   NA         NA         NA        NA  -89.49542  86.02403   19
7   7   52 0.64265299 0.59516900 0.8185208  133.40258  31.93410   17
8   7   70 0.14084638 0.08905376 0.1137092  133.40258  31.93410   17
9   9  120 0.33380777 0.32627085 0.6227753   38.73421 -32.76827   11
10  9  119 0.57976869 0.11692465 0.3070686   38.73421 -32.76827   11
11  9  101 0.76347539 0.94087637 0.4202294   38.73421 -32.76827   11
12  9  102 0.31022532 0.17236697 0.5379526   38.73421 -32.76827   11
13 10   39 0.02508786 0.63041189 0.5683433 -131.46412  35.84172   10

Also in your code, I would clean up where you set your workspace (this should be outside of the loop) and don't use getwd in your loop either, set a variable for the path you need.

  • Thank you for helping me learn! However, in my case I am dealing with hundreds of images that have different extents therefore i can not stack them if I am correct. Furthermore, they are in different directories - this is the reason i am setting the directory within the loop. That was the only way i knew how to import these rasters from different directories, and this also allowed me to take the raster names and put them in columns of the dataframe for the further usage. – davor korman Nov 7 '17 at 16:38
  • Yeah the code above is for the actual process once you identify your rasters in the directory. You can't stack them if they're from different extents but depending on how huge this dataset is perhaps you could stitch the rasters together? Your working directory shouldn't change in a loop, what you want to change is the path variable. I'll edit to clarify. – GISKid Nov 7 '17 at 19:39

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.