1

I've been struggling a lot lately to produce a map in R with the ggplot2 package. So far this is the map I managed to achieve (where each little plygon stands for a city) :

enter image description here

The data used for this map is here, the code used for this map is:

#require 
library(ggplot2) 
library(rgdal)
library(plyr)
library(dplyr)

#open df
data = read.table("datafr_comm.txt", header=T, sep="\t", quote="", dec=".", colClasses=c(rep("factor",)))
data$score <- as.numeric(as.character(data$score))

#open shapefile 
mapa <- readOGR(dsn="France",layer="codes_postaux_region")
mapa <- spTransform(mapa, CRS("+proj=longlat +datum=WGS84"))
#prepare for merging dataframe/shapefile by ID
mapa@data$id <- rownames(mapa@data)
mapa@data   <- join(mapa@data, data, by="ID")
mapa.df     <- fortify(mapa)
mapa.df     <- join(mapa.df,mapa@data, by="id")
#merge dataframe/shapefile
plotDatafr <- left_join(mapa.df, data)

#plot map
ggplot() +
geom_polygon(data = plotDatafr, aes(x = long, y = lat, group = group), fill = "white", colour = NA) +
geom_polygon(data = plotDatafr, aes(x=long, y = lat, group = group, fill=item), colour = NA) + 
scale_fill_manual(values = c("#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7"), name = "") + 
coord_map()

Now, I can't figure out a way to fill automatically the white zones on the map by calculating an average of the surrounding points. I tried many things including the stat_density_2d or stat_bin2d functions, but no way to achieve what I really want.

  • I'm aware of this post but could not adapt it on my own problem... – Mathieu Avanzi Nov 27 '16 at 17:17
  • Is this some sort of discrete classification? So there's no interpolation possible? Do you just want to assign each unknown region to one of the categories? I can't think of a statistically rigorous way of doing it without a complex multivariate spatial model. Plenty of ad-hoc methods spring to mind, like building an adjacency graph and using a colouring scheme based on majority neighbours or something... Either way its not a heatmap and its not smoothing (unless interpolation is possible). – Spacedman Nov 27 '16 at 18:26
  • @spacedman I m more looking for an ad hoc solution since this is kind of a discrete classification (values are assigned to each region). – Mathieu Avanzi Nov 27 '16 at 20:52
  • 1
    So would setting the missing values to the value of the nearest non-missing value be okay? – Spacedman Nov 27 '16 at 21:12
  • If you just want the map to "look nice", you might consider some "photoshop" style display filters instead of trying to get R to generate something? – barrycarter Nov 27 '16 at 22:03
2

Here's the "assign value of nearest non-missing value by centroid distance" solution.

I'll use a simpler data set so you can see what's happening.

require(raster)
require(spdep)
columbus <- shapefile(system.file("etc/shapes/columbus.shp",package="spdep"))

To which I'll add a THING column which is a factor with a lot of missing data. Set the seed for reproducibility:

set.seed(310366);columbus$THING = as.factor(sample(c("A","B","C",NA,NA,NA,NA),49,TRUE))
spplot(columbus,"THING")

original data with missing values

Now I need the FNN package for nearest neighbours. I'll then calculate the nearest neighbours of the centroids of the polygons with missing THING to the ones with non-missing THING:

require(FNN)
srcs = coordinates(columbus)[!is.na(columbus$THING),]
dsts = coordinates(columbus)[is.na(columbus$THING),]
nn = get.knnx(srcs, dsts, 1)

Now nn$nn.index is an index into the non-missing THING values, so I can fill in the missing things with a lookup into the missing things. Like this:

columbus$THING[is.na(columbus$THING)] = columbus$THING[!is.na(columbus$THING)][nn$nn.index]

spplot(columbus,"THING")

nearest-neighbour fill in

To use this correctly, make sure your geospatial data is on a metric coordinate system, so before you project to lat-long. You should probably also wrap it up into a nice function.

fillin <- function(shapes, name){

    srcs = coordinates(shapes)[!is.na(shapes[[name]]),]
    dsts = coordinates(shapes)[is.na(shapes[[name]]),]
    nn = get.knnx(srcs, dsts, 1)

    shapes[[name]][is.na(shapes[[name]])] = shapes[[name]][!is.na(shapes[[name]])][nn$nn.index]
    shapes
}

Then you can just do:

 columbus = fillin(columbus, "THING")

Here's your french data (plotted with spplot rather than faffing around with ggplot, so different colours):

before and after

This was done on mapa just after the data is joined to it, and without doing an spTransform.

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.