# Associating country to grid cell in R?

I’m trying to solve the following problem.
I have a regular grid of points covering Europe

``````nx <- 361 ; ny <- 181
xmin <- -30.0 ; ymin <- 25.0
dx <- 0.25 ; dy <- 0.25
lat <- seq(ymin, ymin+(ny-1)*dy, dy)
lon <- seq(xmin, xmin+(nx-1)*dx, dx)
``````

What I need doing is to associate the country to each grid cell and do country operations.
For example, I have to multiply the population density of France (Spain, Poland, etc) by a given factor.
I have produced an ugly script for masking the grid coordinates with country’s names, but it is not very practical and was wondering if there is a more efficient way to do it.
Here is the code I have produced :

``````m <- matrix(1, ncol=nx, nrow=ny) # create grid with unit (or whatever) value
n <- vector(mode='character', length=nx*ny) # empty character vector with    same length as raster of grid
EUstates <- c('AUT', 'BEL', 'BGR', 'HRV', 'CYP', 'CZE', 'DNK', 'EST', 'FIN', 'FRA', 'DEU', 'GRC',
'HUN', 'IRL', 'ITA', 'LVA', 'LTU', 'LUX', 'MLT', 'NLD', 'POL', 'PRT', 'ROU', 'SVK',
'SVN', 'ESP', 'SWE', 'GBR')
r = raster( m, xmn=xmin,xmx=max(lon),ymn=ymin,ymx=max(lat)) # create a raster of grid 'm' to be used for masking
# now loop through EU countries for masking
for (c in EUstates){
rr <- mask(r, EUc) # assigns NAN to whichever is outsided country c
n[which(!is.na(rr@data@values))] <- c  # assign country name to cells     that are not NAN
}
n[which(n=='')] <-'SEA'  # all remaining empty cells assigned as 'SEA'
``````
• Please do not include thanks and comments of appreciation in your questions - that is what upvotes and the Accept checkmark are for.
– PolyGeo
Feb 16, 2017 at 20:33

Here's a version with a few tweaks:

``````library(sp)
library(raster)
library(rgdal)
library(rmapshaper)
library(tidyverse)

nx <- 361 ; ny <- 181
xmin <- -30.0 ; ymin <- 25.0
dx <- 0.25 ; dy <- 0.25
lat <- seq(ymin, ymin+(ny-1)*dy, dy)
lon <- seq(xmin, xmin+(nx-1)*dx, dx)

# make an empty grid instead so NA = Ocean
m <- matrix(NA, ncol=nx, nrow=ny)
r <- raster(m, xmn=xmin,xmx=max(lon),ymn=ymin,ymx=max(lat))

EUstates <- c('AUT', 'BEL', 'BGR', 'HRV', 'CYP', 'CZE', 'DNK', 'EST',
'FIN', 'FRA', 'DEU', 'GRC', 'HUN', 'IRL', 'ITA', 'LVA',
'LTU', 'LUX', 'MLT', 'NLD', 'POL', 'PRT', 'ROU', 'SVK',
'SVN', 'ESP', 'SWE', 'GBR')

# now loop through EU countries and dump each into a list
for (c in EUstates){
}

# convert list to single spdf
``````

Now simplify the geometry, the linework is unnecessarily detailed for this procedure and its slowing things down:

``````GADM0_Eur_simple <- ms_simplify(GADM0_Europe, keep = 0.05)
``````

Effectively, this next part fills your empty raster with values from column 'ID_0' in the GADM attribute table, which appears to be the numeric ID version of the 'ISO' code. Note that the pixel value is only updated if its center falls on a polygon, so it might be a bit too conservative around coastlines. Easiest solution is to make a raster with smaller cells if that's causing problems. You can always generalise it back to big cells later.

``````GADM0_Europe_raster <- rasterize(GADM0_Eur_simple, r,
field      = 'ID_0',
background = -9999,
update     = TRUE)
``````

Your raster can now be converted to a categorical type and your country name codes can be added to the raster attribute table.

``````GADM0_Europe_raster <- ratify(GADM0_Europe_raster)
rename(ID_0 = ID) %>%
left_join(., GADM0_Eur_simple@data[, c('ID_0', 'ISO')], by = 'ID_0') %>%
# not sure why this got factorised during ms_simplify!
mutate(ISO = as.character(ISO)) %>%
# this last rename isn't optional
rename(ID = ID_0)

You can pull a copy of the RAT out again as a dataframe and do all the calculations you want from there, just be careful copying data back in. Pay attention to `?ratify`. Also, some spatial R functions wipe out the RAT for mystery reasons, which can be frustrating.
• many thanks for the useful hints. Unfortunately I'm unable to install the packages rmapshaper and tidyverse...I got to the point `GADM0_Europe <- do.call('rbind', GADM0)` which gives the message `Error in validObject(res) : invalid class “SpatialPolygons” object: non-unique Polygons ID slot values`