I have a raster file from the FAO's Global Agro-ecological Zone database ( http://gaez.fao.org/Main.html#) which gives the yields of specific crops in 5' grid cells. I would like to assign a country name to each of these grid cell (obviously drop the cells that have no yield data or lie in the oceans). I was thinking of using GADM but I don't know how to reconcile these two sources.

Ultimately, I want to create a data frame with the following columns:

Grid cell ID | Longitude | Latitude | Yield (from raster) | Country

I'm quite clueless as to where to even begin. For some broader context, I'm trying to follow this paper: https://economics.mit.edu/files/7536 by summing up the crop production in each country from the grid cells.

closed as off-topic by ahmadhanb, LaughU, Jochen Schwarze, nmtoken, whyzar Oct 31 '18 at 13:11

This question appears to be off-topic. The users who voted to close gave this specific reason:

  • "When seeking help to debug/write/improve code always provide the desired behavior, a specific problem/error and the shortest code (as formatted text, not pictures) needed to reproduce it in the question body. Providing a clear problem statement and a code attempt helps others to help you." – ahmadhanb, LaughU, Jochen Schwarze, nmtoken, whyzar
If this question can be reworded to fit the rules in the help center, please edit the question.


See ?raster::extract or ?raster::crosstab. For example:

# example regions
p <- shapefile(system.file("external/lux.shp", package="raster"))
# example yield raster
yield <- raster(p)
dim(yield) <- c(100, 100)
values(yield) <- 1:ncell(yield)
yield[sample(ncell(yield))[1:1000]] <- NA

To get yield by polygon (e.g. country), you can use extract and combine with the region names

e <- extract(r, p, fun=mean, na.rm=T)
data.frame(name=p$NAME_2, yield=e)
#               name    yield
#1          Clervaux 1493.427
#2          Diekirch 4307.677
#3           Redange 5168.996
#4           Vianden 3446.380

To get yield by grid cell, you can do

id <- rasterize(p, yield)
yield <- mask(yield, id)
s <- stack(id, y)
a <- rasterToPoints(s)
cells <- cellFromXY(s, a[,1:2])
d <- data.frame(ID=cells, lon=a[,1], lat=a[,2], yield=a[,4], country=p$NAME_2[a[,3]])
#   ID      lon      lat yield  country
#1  35 6.014659 50.17795    35 Clervaux
#2  36 6.022500 50.17795    36 Clervaux
#3 129 5.967612 50.17061   129 Clervaux
#4 130 5.975453 50.17061   130 Clervaux
#5 131 5.983294 50.17061   131 Clervaux
#6 132 5.991135 50.17061   132 Clervaux
  • Thank you so much! This worked beautifully even though I didn't really provide you with any of my dataset. – nicholasflamel Nov 2 '18 at 2:20

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