How to use function crop(x, y, ...) from raster package so that the output is information about y (raster) contained in x (shapefile, extent object)

Question

I am kind of new to R and especially to working with GIS in R, so I hope I explain this right.

I want to get the function crop (x,y...) from package raster to merge/overlay (not sure what's the correct expression to use) a raster file with a shapefile. Essentially, the shapefile is a 5 x 5 km grid of Sweden and the raster is a landuse raster of Sweden. By using the crop function I want to get a result which, for each 5x5 km square from the shapefile, gives information extracted from the raster file about the type of landuse in each square.

However, when I use crop and then write.csv from what I cropped 'crop (landuse, ekrut)" (see code below), I basically just get the ekrut (shapefile) as output csv file, there are no columns added of the landuse raster indicating which square has which land use class.

I am sorry but I cannot share the actual shapefile here, if it will make a difference somehow, the landuse data can be downloaded from here: http://gpt.vic-metria.nu/data/land/NMD/NMD2018_basskikt_ogeneraliserad_Sverige_v1_0.zip

Here's what I tried to do so far

library (raster)
library (rgdal)
library (sp)

this is the coordinate system/projection for the GIS file. Each coordinate system has an epsg code (http://spatialreference.org/).

sweref.def <- "+init=epsg:3006" 

# read in shapefile
ekrut <- readOGR (dsn   = "//storage-al.slu.se/student$/nilc0001/Desktop/Nina/Ekrut", 
                  layer = "ekrut_5x5_flat", 
                  p4s   = sweref.def) 

ekrut
# class       : SpatialPolygonsDataFrame 
# features    : 19192 
# extent      : 266333, 921700, 6131565, 7675329  (xmin, xmax, ymin, ymax)
# coord. ref. : +init=epsg:3006 +proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
# variables   : 7
# names       :    AREA, PERIMETER,  GGD_, GGD_ID,     BK, Bk_num, BK_flat 
# min values  : 2.5e+07,     20000,     2,      0, 10A 0f, 012 77,   10A0f 
# max values  : 2.5e+07,     20000, 19208,  19230,  9J 9d, 329 48,    9J9d 

#read in raster
landuse <- raster ("nmd2018bas_ogeneraliserad_v1_0.tif")

landuse

# class       : RasterLayer 
# dimensions  : 157992, 71273, 11260563816  (nrow, ncol, ncell)
# resolution  : 10, 10  (x, y)
# extent      : 208450, 921180, 6091140, 7671060  (xmin, xmax, ymin, ymax)
# coord. ref. : +proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
# data source : //storage- al.slu.se/student$/nilc0001/Desktop/Nina/landuse/nmd2018bas_ogeneraliserad_v1_0.tif 
# names       : nmd2018bas_ogeneraliserad_v1_0 
# values      : 0, 128  (min, max)
# attributes  :
#         ID      COUNT Opacity Klass
#  from:   0 5204803484       0      
#  to  : 255          0       0   

#first few rows of attribute table of landuse
levels (landuse)   

# [[1]]
#      ID      COUNT Opacity                                              Klass
# 1     0 5204803484       0                                                   
# 2     1          0     255                                                   
# 3     2  382320369     255                                      Öppen våtmark
# 4     3  258249590     255                                           Åkermark

# crop and write csv
landuse.ekrut <- crop (landuse, ekrut)
write.csv (landuse.ekrut,"landuse.ek.csv")

Like I said, when I use crop and then write.csv from what I cropped I basically just get the ekrut (shapefile) as output csv file, there are no columns added of the landuse raster indicating which square has which land use class.

I would like to have a .csv which, for each square (there are 19 191 of them), gives me information on what type of landuse is present in this square.

If there is a better way to approach this, please indicate. I hope I was clear with my explanation!

Thanks.

Answer 1

Even if your code was working, still you could not export this output into CSV unless you use something like majority filter. This is basically because of the fact that one polygon may contain more than one land-use.

One way to do this is as below:

library(raster)
library(stats)

#set.seed to have a reproducible example
set.seed(1523)

#making a raster layer
r <- raster(ncol=36, nrow=18, vals=floor(runif(648, 5.0, 7.5)))
rr <- aggregate(r, fact=3)
f.net <- rasterToPolygons(rr) # transform into a polygon

#simple ovelaying plot for visualisation
plot(r, main="Raster and the overlaying poly")
plot(f.net, col=rgb(1, 0, 0, alpha = 0.3), add=TRUE)

ext.r <- extract(r, f.net) # returns values per ploy

# a simple mode filter
mode <- function(x){ 
  ta = table(x)
  tam = max(ta)
  mod = as.numeric(names(ta)[ta == tam])
  mod = mod[1] 
  return(mod)
}

# each poly contains several values
head(ext.r)
# [[1]]
# [1] 6 5 5 6 5 6 5 6 5
# 
# [[2]]
# [1] 5 5 6 7 6 7 6 5 7
# 
# [[3]]
# [1] 5 5 7 6 6 7 6 6 7
# 
# [[4]]
# [1] 6 7 5 5 7 6 5 6 5
# 
# [[5]]
# [1] 5 6 5 5 5 5 6 5 5
# 
# [[6]]
# [1] 7 5 6 6 6 5 5 6 5


# a mode filter can be applied to pas majority values to each poly
mode.ext.r <- lapply(ext.r, mode)
head(mode.ext.r)
# [[1]]
# [1] 5
# 
# [[2]]
# [1] 5
# 
# [[3]]
# [1] 6
# 
# [[4]]
# [1] 5
# 
# [[5]]
# [1] 5
# 
# [[6]]
# [1] 5
write.table(mode.ext.r, file = 'mode_ext_r.csv', row.names = FALSE, na = '')

Answer 2

You can extract all classes for each polygon, make a table for each polygon, and then combine them.

Building off on Majid's example data.

library(raster)
set.seed(1523)
r <- raster(ncol=9, nrow=6, vals=sample(10, 9*6, replace=TRUE))
rr <- aggregate(raster(r), fact=3)
pols <- rasterToPolygons(rr) # transform into a polygon
plot(r); lines(pols);  text(pols, 1:length(pols))

Extract, create table, and combine

e <- extract(r, pols)
x <- lapply(e, table)
y <- lapply(1:length(x), function(i) data.frame(pol=i, as.data.frame(x[[i]])))
z <- do.call(rbind, y)

head(z, 10)
#   pol Var1 Freq
#1    1    1    2
#2    1    2    2
#3    1    3    2
#4    1    6    2
#5    1    9    1
#6    2    1    1
#7    2    2    1
#8    2    4    2
#9    2    6    2
#10   2    7    1

And save to csv

write.csv(z, "test.csv", row.names=FALSE)