[adegenet-forum] interpolation in a map

[Jombart, Thibaut]

Zhongkui.Luo at csiro.au wrote:
> Hi,
>  
> Sorry for the bother. A R problem looks for your suggestion.
>  
> I want to draw a figure with filled contour in an irregular map.
>  
> But I found that the contours overlap the boundaries of the map, and also some regions in the map were not filled.
>  
> I used the following R code:
>  
>    
> Map<-readShapeSpatial("data1") # data1, a shapefile, the whole Australian map with boundary.
> plot(Map)
> data2.li<-interp(data2$x,data2$y,data2$z) #data2, txt file with three columns: x, longitude; y, latitude; z, variable (e.g., rainfall)
> filled.contour(data2.li, color=topo.colors, add=TRUE)
>  
>  
> Are there any R functions can be used to limit the filled contours just inside the map (no overlap and no blank)?
>  
> Do you have any ideas to deal with this problem? Thanks very much for that!
>  
> Best,
>  
>  
>
> Zhongkui Luo
>
> PhD candidate
>
> CSIRO Land & Water, Butler Lab.
>
> Clunies Ross Street, Acton
>
> Canberra, ACT 2601, Australia
>
> Ph: +61 2 6246 5594  E_Mail: Zhongkui.luo at csiro.au
>
>  
>

Dear Luo,

I am afraid your post is off-topic for this ML. Also note that it is generally considered as impolite to post the same question at the same time to different mailing lists, as you did here. Lastly, you could find your answer in one command line:
RSiteSearch("blank map contour")

which takes you to relevant results found in R mailing lists archives, where the solution has been posted by Duncan Murdoch last year.

However, that kind of issue may occur to people analysing genetic data, for instance when geographically mapping the results of a multivariate analysis. Here is a ready-to-use solution:

##########
library(akima)
library(maps)
library(mapdata)

GRID.RES <- 500

## original values
x <- seq(100,170,le=10) + runif(10,-5,5)
y <- seq(-50,0,le=10) + runif(10,-5,5)
z <- x*y + rnorm(length(x),2)

## interpolation using akima
x0 <- seq(100,170,le=GRID.RES)
y0 <- seq(-50,0,le=GRID.RES)
z.hat <- interp(x,y, z, xo=x0, yo=y0, linear = FALSE, extrap=TRUE)

## find loc in australia, set others to NA
myGrid <- expand.grid(x0, y0)
temp <- map.where(x=myGrid[,1], y=myGrid[,2])
toKeep <- grep("australia",temp, ignore.case=TRUE)
toRemove <- setdiff(1:length(z.hat$z), toKeep)
z.hat$z[toRemove] <- NA

## plot
image(z.hat)
map(add=TRUE, lwd=3)
contour(z.hat,add=TRUE)
############

The quality of the result depends on the grid size (GRID.RES); GRID.RES=1000 takes a few seconds on my computer but provides publication-quality figures. For genetic data analysis, 'z' could be values of a principal component to be mapped, while x and y would be spatial coordinates of samples.

Best,

Thibaut.

-- 
######################################
Dr Thibaut JOMBART
MRC Centre for Outbreak Analysis and Modelling
Department of Infectious Disease Epidemiology
Imperial College - Faculty of Medicine
St Mary’s Campus
Norfolk Place
London W2 1PG
United Kingdom
Tel. : 0044 (0)20 7594 3658
t.jombart at imperial.ac.uk
http://www1.imperial.ac.uk/medicine/people/t.jombart/
http://adegenet.r-forge.r-project.org/