4

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){
    EUc <- getData("GADM", country=c, level=0)
    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 '17 at 20:33
4

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
GADM0 <- list()
for (c in EUstates){
  GADM0[[c]] <- getData("GADM", country=c, level=0)
}

# convert list to single spdf
GADM0_Europe <- do.call('rbind', GADM0)

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)
rat <- levels(GADM0_Europe_raster)[[1]] %>%
  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)

levels(GADM0_Europe_raster)[[1]] <- rat

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 – efz Feb 16 '17 at 16:18
  • I'm pretty sure the error you describe was patched quite some time ago - its discussed here gis.stackexchange.com/questions/32732/…. Is a firewall the reason you can't update or install anything? – obrl_soil Feb 16 '17 at 22:24

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.