[R-sig-Geo] FW: interpolation in a map

Wed Sep 23 02:08:19 CEST 2009

Dear all,

Thanks very much for Thibaut's help first!

Dr. Thibaut and Other colleagues have gave a well solution for drawing filled contours in an irregular map (See the following emails).
However, the drawing region looks like must be already in the map database, such as the whole USA  and AUS. 
If the region is not in the initial map, how to use this region's shapefile to get the map with filled contours?

The data might be like this:

Data1: a shapefile call "xxx.shp", e.g., a catchment.
Data2: a txt file with three columns: x, longitude; y,latitude;z, variable (e.g., rainfall).

This work would like be very easy in GIS, but how is it done in R?

Thanks very much for that!

Best,

Zachary

-----Original Message-----
From: Jombart, Thibaut [mailto:t.jombart at imperial.ac.uk] 
Sent: Tuesday, 22 September 2009 11:21 PM
To: Luo, Zhongkui (CLW, Black Mountain)
Subject: RE: interpolation in a map

Hi Zachary (or Zhongkui?), 

I think I would make the proper shapefile in a GIS, then import it in R using maptool and use the same method as for Australia. It it also possible, I guess, to build your polygon by hand in R, but in think much more cumbersome.

Note that a ML is specially dedicated to this kind of issue:
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Best,

Thibaut.
________________________________________
From: prvs=5094c3904=Zhongkui.Luo at csiro.au [prvs=5094c3904=Zhongkui.Luo at csiro.au] On Behalf Of Zhongkui.Luo at csiro.au [Zhongkui.Luo at csiro.au]
Sent: 22 September 2009 10:53
To: Jombart, Thibaut
Subject: RE: interpolation in a map

Dear Thibaut,

Sorry to bother you again!

Your method can work well at whole Australia. That's just I want.
But how can I conduct the similar process at a specific region inside Australia.
It looks like that Australia is not divided into different subregions, such as states or counties, and just a whole map of AUS in the package.
So there are no other locs.

At this situation, the question is:
How to add subregions to a map according to corresponding shapefiles?
Does R has specific function to do this? Or do you have some ideas or information to deal with it?

Thanks very much for that.

Best,

Zachary

-----Original Message-----
From: Jombart, Thibaut [mailto:t.jombart at imperial.ac.uk]
Sent: Sunday, 20 September 2009 10:54 PM
To: Luo, Zhongkui (CLW, Black Mountain); adegenet-forum at lists.r-forge.r-project.org
Subject: RE: interpolation in a map

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/