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'
PolyGeo
  • 65,136
  • 29
  • 109
  • 338
efz
  • 143
  • 4
  • 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

1 Answers1

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.

obrl_soil
  • 3,752
  • 15
  • 34
  • 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 https://gis.stackexchange.com/questions/32732/proper-way-to-rbind-spatialpolygonsdataframes-with-identical-polygon-ids. Is a firewall the reason you can't update or install anything? – obrl_soil Feb 16 '17 at 22:24