Wednesday 17 February 2016

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){
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'

Answer




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.


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...