7

I want to create a contour of variable z with the x,y,z data. However, it seems like we need to provide the data in increasing order.

I tried to use some code but it gave me the error.

I tried the following code: Trial 1:

age2100 <- read.table("temp.csv",header=TRUE,sep=",")

x <- age2100$x
y <- age2100$y
z <- age2100$z

contour(x,y,z,add=TRUE,col="black")

I got the following error

Error in contour.default(x, y, z, add = TRUE, col = "black") : increasing 'x' and 'y' values expected

I then tried to use ggplot2 to create the contour. I used the following code:

library("ggplot2")
library("MASS")
library("rgdal")
library("gpclib")
library("maptools")
age2100 <- read.table("temp.csv",header=TRUE,sep=",")
v <- ggplot(age2100, aes(age2100$x, age2100$y,z=age2100$z))+geom_contour()
v

I got the following error:

Warning message:

Not possible to generate contour data

Please find the data on the following location https://www.dropbox.com/s/mg2bo4rcr6n3dks/temp.csv

Can anybody tell me how to create the contour data from the third variable (z) from the temp.csv ? I need to do these many times so I am trying to do on R instead of Arcgis.

agstudy
  • 119,832
  • 17
  • 199
  • 261
Jd Baba
  • 5,948
  • 18
  • 62
  • 96
  • Just a long shot, have you tried `age2100 <- age2100[with(age2100, order(x, y)),]` – sebastian-c Jan 17 '13 at 09:07
  • @sebantian-c I just used your suggestion. I checked the data and it is sorted but still I am not able to get the graph. I then tried both methods described above. I still get the same error. – Jd Baba Jan 17 '13 at 09:21
  • I believe you need a full matrix - you can interpolate one with your xyz data using the `interp`function of the package `akima`. – Marc in the box Jan 17 '13 at 09:49
  • 1
    @Jdbaba for me it is a ggplot2 bug. for some reason stat_contour can't generate contourLines. – agstudy Jan 17 '13 at 10:11

1 Answers1

9

Here is an example of how one interpolates using interp from the akimapackage:

age2100 <- read.table("temp.csv",header=TRUE,sep=",")

x <- age2100$x
y <- age2100$y
z <- age2100$z

require(akima)

fld <- interp(x,y,z)

par(mar=c(5,5,1,1))
filled.contour(fld)

enter image description here

Here is an alternate plot using the imagefunction (this allows some flexibility to adding lower level plotting functions (requires the image.scale function, found here):

source("image.scale.R") # http://menugget.blogspot.de/2011/08/adding-scale-to-image-plot.html

x11(width=5, height=6)
layout(matrix(c(1,2), nrow=1, ncol=2), widths=c(4,1), height=6, respect=TRUE)
layout.show(2)

par(mar=c(4,4,1,1))
image(fld)
contour(fld, add=TRUE)
points(age2100$x,age2100$y, pch=".", cex=2)

par(mar=c(4,0,1,4))
image.scale(fld$z, xlab="", ylab="", xaxt="n", yaxt="n", horiz=FALSE)
box()
axis(4)
mtext("text", side=4, line=2.5)

enter image description here

Marc in the box
  • 11,769
  • 4
  • 47
  • 97
  • Marc in the box, Thank you so much for your solution. The plot looks really nice. I have one question. When I add the points on the above plot it seems the points are little bit off the contour lines that is created. I may be wrong. Is it possible to add the points in the above plot ? I used the following code after the plot points(age2100$x,age2100$y,pch=0). Is it possible to show the lines also in the plot ? Just like assigning the contour interval and then show will lines. – Jd Baba Jan 17 '13 at 10:53
  • Yes, the lattice function `filled.contour`(and others) split the device region, and thus you cannot directly add lower level plots to the plot only. It's not pretty, but I added a few lines to help you with these additions. I constantly use my own function for adding the color scale for precisely this reason that you mention. – Marc in the box Jan 17 '13 at 11:19
  • Thank you so much for editing the code. I was just wondering whether it is possible to add a shapefile and then clip the part of contour data and then only show the contour within the shapefile ? Can you give me some suggestions on how that can be achieved ? – Jd Baba Jan 18 '13 at 02:27
  • Marc in the box, All the solution in perfect but I am confused about the legend. In the first part, the legend is good but when you modified the code with the points and contour added , then the legend has changed. Can you give suggestions on how to change to the right scale ? – Jd Baba Jan 18 '13 at 03:03
  • Hi @Jdbaba - I fixed the legend units. It just needed the input of `fld$z` instead of `fld` to get the range correct. I am unsure how to clip the contour plot. You could try overlaying some other plot to shade out the region that you don't want to see (e.g. `image(, add=TRUE)`) – Marc in the box Jan 18 '13 at 12:57