Extracting RAMAS GIS carrying capacities (.KCH) with R

RAMAS GIS has some really powerful capabilities with PVA, but a lot of the model output is hidden amongst many files. One thing that annoys me is the default map in RAMAS. It shows the carrying capacity of each population (presumably at the beginning of the simulation?), without any spatial context.

I wasn't happy with this, so I decided to extract the data and create a map of my own. 

The carrying capacity information is held in two locations. Stable carrying capacities are held within the main .MP file. Metapopulations with changing carrying capacities store the carrying capacity of each year in a .KCH file, which the .MP file refers to. There are also linearly changing carrying capacities (gradient kept in the .MP file), but I've ignored them for now.

This code will extract carrying capacities at the end of the simulation from the .MP and .KCH files (kept in a folder designated by scenario). It then plots populations onto a world map.

library(stringr)
library(maps)

rm(list = ls())
setwd("C:/Evan/R/ramas/kchextract") #Working directory should hold folders with each PVA scenario
scenario <- "BFS_IP85_NA/" # Folder holding .MP and .KCH files for one scenario
files <- list.files(scenario)

mp_file <- read.table(paste(scenario, subset(files, str_sub(files, -3, -1) == ".MP"), sep=""), header = FALSE, sep = "@")

mp_data <- colsplit(mp_file[33 : (which(mp_file == "Migration")-1), ], split = ",", names = c(1:27))
mp_data$X10 <- as.character(mp_data$X10)
mp_data$X10 <- ifelse(mp_data$X10 == "" | is.na(mp_data$X10), mp_data$X7, mp_data$X10)

get_k <- function(x) {
  if(str_sub(mp_data$X10[x], -3, -1) == "KCH") {
    read.table(paste(scenario, mp_data$X10[x], sep = ""))[nrow(read.table(paste(scenario, mp_data$X10[x], sep = ""))),]
  } else {
    as.numeric(mp_data$X7[x])
  }
}

mp_data$k <- aaply(c(1:length(mp_data$X10)), 1, get_k)

## This converts the pixel locations to true long/lats. Find corner of your map to convert yourself.
mp_data$x.long <- (mp_data$X2 / 7200 * 60) + 90 
mp_data$y.lat <- 53.5619398133 - (mp_data$X3 / 5877 * 48.97500168)

world <- map_data('world')

p <- ggplot(world, aes(x = long, y = lat)) +
  geom_polygon(aes(group = group), fill = "grey98", alpha = 0.5, colour = "black") +
  geom_point(data = mp_data, aes(x = x.long, y = y.lat, size = k), alpha = 0.3) +
  scale_size_continuous(range = c(1,12)) +
  guides(size=F, alpha=F) +
  coord_cartesian(xlim = c(95, 145), ylim = c(4.5869381328, 40)) +
  xlab("Longitude") +
  ylab("Latitude") +
  theme_evp() +
  theme(legend.justification=c(0, 1), legend.position=c(0, 1)) +
  scale_alpha_continuous(range=c(0, 1))


ggsave(filename = paste(str_sub(scenario, 1, -2), ".jpeg", sep = ""), plot = p, width = 35.88, height = 23.8, units = "cm")

Note: I realised whilst writing this blog post that there is a package called mptools that will do some of this stuff. Don't think it replicates this exactly, but oh well. Hopefully this is still useful to some.