Re: [R-sig-Geo] 3rd+ rasters added using plot(r, add=TRUE) are misregistered
Good catch. I have no solution -- but I get the same problem on an earlier version of R in a different environment (so it is not just because of your setup). My session info is below. Cheers, Gareth. sessionInfo() R version 3.0.1 (2013-05-16) Platform: x86_64-unknown-linux-gnu (64-bit) locale: [1] LC_CTYPE=en_AU.UTF-8 LC_NUMERIC=C [3] LC_TIME=en_AU.UTF-8LC_COLLATE=en_AU.UTF-8 [5] LC_MONETARY=en_AU.UTF-8LC_MESSAGES=en_AU.UTF-8 [7] LC_PAPER=C LC_NAME=C [9] LC_ADDRESS=C LC_TELEPHONE=C [11] LC_MEASUREMENT=en_AU.UTF-8 LC_IDENTIFICATION=C attached base packages: [1] grid stats graphics grDevices utils datasets methods [8] base other attached packages: [1] rgdal_0.8-10rgeos_0.2-19raster_2.1-59 maptools_0.8-25 [5] lattice_0.20-15 sp_1.0-11 foreign_0.8-53 -- View this message in context: http://r-sig-geo.2731867.n2.nabble.com/3rd-rasters-added-using-plot-r-add-TRUE-are-misregistered-tp7585262p7585271.html Sent from the R-sig-geo mailing list archive at Nabble.com. ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] covariance matrix
Dear Erin, you might want to take a look at the gstat function variogramLine and set the parameter covariance to TRUE. The argument dist_vector takes as well a distance matrix and will then return a covariance matrix. HTH, best wishes, Ben Hodgess, Erin hodge...@uhd.edu wrote: Hi! Here is a goofy question, please: how do I get the covariance matrix based on the estimated model, say exponential, please? I was looking at the variogram function in gstat and I think that the gamma values may be part of the solution, but I'm not sure how to finish it out. I'm experimenting with putting together a toy data set through the variogram and kriging process by hand (to a certain extent) Thanks for any help! Sincerely, Erin [[alternative HTML version deleted]] ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo -- Benedikt Gräler Institute for Geoinformatics University of Münster http://www.ifgi.de/graeler +49 251 83 33082 ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] covariance matrix
Erin, http://ifgi.uni-muenster.de/~epebe_01/mstp/lec5.html http://ifgi.uni-muenster.de/~epebe_01/mstp/lec7.html contain some minimal, one-dimensional kriging functions that I used in my class recently. R markdown source files are on https://github.com/edzer/mstp On 12/05/2013 04:36 AM, Hodgess, Erin wrote: Hi! Here is a goofy question, please: how do I get the covariance matrix based on the estimated model, say exponential, please? I was looking at the variogram function in gstat and I think that the gamma values may be part of the solution, but I'm not sure how to finish it out. I'm experimenting with putting together a toy data set through the variogram and kriging process by hand (to a certain extent) Thanks for any help! Sincerely, Erin [[alternative HTML version deleted]] ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo -- Edzer Pebesma Institute for Geoinformatics (ifgi), University of Münster Heisenbergstraße 2, 48149 Münster, Germany. Phone: +49 251 83 33081 http://ifgi.uni-muenster.de GPG key ID 0xAC227795 signature.asc Description: OpenPGP digital signature ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
[R-sig-Geo] Strange lines when projecting a world map
Dear list, as in the title, I found that when projecting a world map by mollweide projection and setting a longitude value, some strange parallels are drawn. The real problem is that these lines âinteractâ with a raster when performing the âoverâ function. How can I manage it? Here I provide a simulated example to let you understand what I mean Thanks. #Example of a world map mollweide projection with a specific longitude value library(sp) library(maptools) library(rgdal) data(wrld_simpl) w_pol-wrld_simpl w_pol-spTransform(w_pol, CRS(+proj=moll +lon_0=66)) plot(w_pol)  # Intersecting a grid with the projected world map gr-expand.grid(x = seq(-180, 180, length.out = 100), y = seq(-90, 90, length.out = 100)) coordinates(gr)-~x+y proj4string(gr)-CRS(+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs) gr-spTransform(gr, CRS(+proj=moll +lon_0=66)) plot(gr[w_pol,])   Dr. Francesco Carotenuto (Ph.D.) Department of Earth Sciences Federico II University Largo S. Marcellino 10, 80138 NAPLES (ITALY) http://francesco-carotenuto.blogspot.it/ http://3dpaleontology.blogspot.it/ [[alternative HTML version deleted]] ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] original coordinates of the redwood data?
On Wed, Dec 4, 2013 at 8:01 PM, Marcelino de la Cruz marcelino.delac...@upm.es wrote: .. The plot of the data published by Strauss (1975) was scanned and digitised by Sandra Pereira, University of Western Australia, 2002. So you can digitize yourself the original true coordinates from the plot in Strauss (1975). Except the plot in Strauss (1975) doesn't have any coordinate labels! Nor a scale bar, north arrow or any of that chart junk. We have a clue to the scale in that the distance 'r1' marked on the plot is 'about six feet'. If anyone would like to rescale the redwood data into units of 'about six feet' then I think that would be a nice waste of their time. Barry ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] covariance matrix
Thanks to all for the good help...it's great! Sincerely, Erin From: r-sig-geo-boun...@r-project.org [r-sig-geo-boun...@r-project.org] on behalf of Edzer Pebesma [edzer.pebe...@uni-muenster.de] Sent: Thursday, December 05, 2013 2:56 AM To: r-sig-geo@r-project.org Subject: Re: [R-sig-Geo] covariance matrix Erin, http://ifgi.uni-muenster.de/~epebe_01/mstp/lec5.html http://ifgi.uni-muenster.de/~epebe_01/mstp/lec7.html contain some minimal, one-dimensional kriging functions that I used in my class recently. R markdown source files are on https://github.com/edzer/mstp On 12/05/2013 04:36 AM, Hodgess, Erin wrote: Hi! Here is a goofy question, please: how do I get the covariance matrix based on the estimated model, say exponential, please? I was looking at the variogram function in gstat and I think that the gamma values may be part of the solution, but I'm not sure how to finish it out. I'm experimenting with putting together a toy data set through the variogram and kriging process by hand (to a certain extent) Thanks for any help! Sincerely, Erin [[alternative HTML version deleted]] ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo -- Edzer Pebesma Institute for Geoinformatics (ifgi), University of Münster Heisenbergstraße 2, 48149 Münster, Germany. Phone: +49 251 83 33081 http://ifgi.uni-muenster.de GPG key ID 0xAC227795 ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Strange lines when projecting a world map
On Thu, 5 Dec 2013, Francesco Carotenuto wrote: Dear list, as in the title, I found that when projecting a world map by mollweide projection and setting a?? longitude value, some strange parallels are drawn. The real problem is that these lines ???interact??? with a raster when performing the ???over??? function. How can I manage it? The underlying question is what you want to do with the objects after projection. The problem in your example is that some polygons/lines cross -180+66 degrees, so get reflected back across. A possible solution is to clip the polygons before projection at -180+66 degrees +/- a small number, but before looking at this in more detail, I need to know what the projected object is intended for. When +lon_0=0, there is no such effect, as this data set is clipped at -180 already. Roger Here I provide a simulated example to let you understand what I mean Thanks.?? #Example of a world map mollweide projection with a specific longitude value library(sp) library(maptools) library(rgdal) data(wrld_simpl) w_pol-wrld_simpl w_pol-spTransform(w_pol, CRS(+proj=moll +lon_0=66)) plot(w_pol) ?? # Intersecting a grid with the projected world map gr-expand.grid(x = seq(-180, 180, length.out = 100), y = seq(-90, 90, length.out = 100)) coordinates(gr)-~x+y proj4string(gr)-CRS(+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs) gr-spTransform(gr, CRS(+proj=moll +lon_0=66)) plot(gr[w_pol,]) ?? ?? Dr. Francesco Carotenuto (Ph.D.) Department of Earth Sciences Federico II University Largo S. Marcellino 10, 80138 NAPLES (ITALY) http://francesco-carotenuto.blogspot.it/ http://3dpaleontology.blogspot.it/?? [[alternative HTML version deleted]] -- Roger Bivand Department of Economics, Norwegian School of Economics, Helleveien 30, N-5045 Bergen, Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43 e-mail: roger.biv...@nhh.no ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Strange lines when projecting a world map
Thank you for the answer. The object is the prediction location in an interpolation process. The function of overlay is involved in an automatic process and the longitude (and the kind of projection too) can change according to the dataset that is, in some cases, worldwide distributed. Dr. Francesco Carotenuto (Ph.D.) Department of Earth Sciences Federico II University Largo S. Marcellino 10, 80138 NAPLES (ITALY) http://francesco-carotenuto.blogspot.it/ http://3dpaleontology.blogspot.it/ Il Giovedì 5 Dicembre 2013 21:23, Roger Bivand roger.biv...@nhh.no ha scritto: On Thu, 5 Dec 2013, Francesco Carotenuto wrote: Dear list, as in the title, I found that when projecting a world map by mollweide projection and setting a?? longitude value, some strange parallels are drawn. The real problem is that these lines ???interact??? with a raster when performing the ???over??? function. How can I manage it? The underlying question is what you want to do with the objects after projection. The problem in your example is that some polygons/lines cross -180+66 degrees, so get reflected back across. A possible solution is to clip the polygons before projection at -180+66 degrees +/- a small number, but before looking at this in more detail, I need to know what the projected object is intended for. When +lon_0=0, there is no such effect, as this data set is clipped at -180 already. Roger Here I provide a simulated example to let you understand what I mean Thanks.?? #Example of a world map mollweide projection with a specific longitude value library(sp) library(maptools) library(rgdal) data(wrld_simpl) w_pol-wrld_simpl w_pol-spTransform(w_pol, CRS(+proj=moll +lon_0=66)) plot(w_pol) ?? # Intersecting a grid with the projected world map gr-expand.grid(x = seq(-180, 180, length.out = 100), y = seq(-90, 90, length.out = 100)) coordinates(gr)-~x+y proj4string(gr)-CRS(+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs) gr-spTransform(gr, CRS(+proj=moll +lon_0=66)) plot(gr[w_pol,]) ?? ?? Dr. Francesco Carotenuto (Ph.D.) Department of Earth Sciences Federico II University Largo S. Marcellino 10, 80138 NAPLES (ITALY) http://francesco-carotenuto.blogspot.it/ http://3dpaleontology.blogspot.it/?? [[alternative HTML version deleted]] -- Roger Bivand Department of Economics, Norwegian School of Economics, Helleveien 30, N-5045 Bergen, Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43 e-mail: roger.biv...@nhh.no [[alternative HTML version deleted]] ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] 3rd+ rasters added using plot(r, add=TRUE) are misregistered
I just wonder if this is related to an existing resize problem when overplotting? I see this on the windows() graphics device, you need to try these two sets of plot commands and resize the window (especially stretch it in one direction) to see the problem at ## 1 and ## 2 v - list(x = seq(0, 1, length = nrow(volcano)), y = seq(0, 1, length = ncol(volcano)), z = volcano) cl - contourLines(v, levels = 160)[[2]] image(v, xaxs = i) lines(cl) ## (1) now resize the window, it's all good - the contour stays in the ## right place ## however, resize after adding the lines to a raster plot ## and the relationship is broken library(raster) plot(raster(v), xaxs = i) lines(cl) ## (2) May be totally unrelated but I thought it worth pointing out. It's hard to trace through the setting up that raster's plot does (well, I find it hard to do - I've tried a few times to figure out where these problems are happening). Cheers, Mike. On Wed, Dec 4, 2013 at 5:57 AM, Josh O'Brien joshmobr...@gmail.com wrote: Hi all, Briefly, running the following code will not produce a single grey ring, as would be expected. The first three raster layers all stack up neatly, but the 4th and subsequent raster layers are all shifted/skewed to the right. It appears that the plotting engine loses track of the fact that some space in the original plot was allocated to a legend, and just begins plotting as if a wider canvas were available to it. I'm not sure whether this is ultimately an issue with raster or sp or something else, which is why I'm noting it here on R-sig-geo. (I'm also fully aware that there are other and better ways to add several raster layers to a plot, but still think the following behavior qualifies as a bug that I ought to report.) library(maptools) ## Only needs to be installed for example data library(raster) library(rgeos) ## Create an example raster p - shapefile(system.file(shapes/co37_d90.shp, package=maptools)) p - p[31,] ## A tall narrow county polygon pr - gDifference(gBuffer(p, width=.01), p) r - rasterize(pr, raster(extent(pr), ncol=100, nrow=100)) ## These three are properly registered on one another plot(r, col=yellow, legend=FALSE) plot(r, col=green, legend=FALSE, add=TRUE) plot(r, col=grey, legend=FALSE, add=TRUE) ## All subsequent layers are improperly shifted/skewed to right plot(r, col=yellow, legend=FALSE, add=TRUE) plot(r, col=blue, legend=FALSE, add=TRUE) plot(r, col=red, legend=FALSE, add=TRUE) plot(r, col=grey20, legend=FALSE, add=TRUE) ## Following the above, SpatialPolygons are also shifted/skewed plot(p) My session info, in case it matters, is: sessionInfo() R version 3.0.2 (2013-09-25) Platform: i386-w64-mingw32/i386 (32-bit) locale: [1] LC_COLLATE=English_United States.1252 [2] LC_CTYPE=English_United States.1252 [3] LC_MONETARY=English_United States.1252 [4] LC_NUMERIC=C [5] LC_TIME=English_United States.1252 attached base packages: [1] stats graphics grDevices utils datasets methods base other attached packages: [1] rgdal_0.8-14 rgeos_0.3-2 raster_2.1-66 sp_1.0-14 loaded via a namespace (and not attached): [1] compiler_3.0.2 grid_3.0.2 lattice_0.20-24 tools_3.0.2 [[alternative HTML version deleted]] ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo -- Michael Sumner Hobart, Australia e-mail: mdsum...@gmail.com ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Strange lines when projecting a world map
On Thu, 5 Dec 2013, Francesco Carotenuto wrote: Thank you for the answer. The object is the prediction location in an interpolation process. The function of overlay is involved in an automatic process and the longitude (and the kind of projection too) can change according to the dataset that is, in some cases, worldwide distributed. This doesn't help much - do you mean that the grid is of prediction locations? Are the polygons used for anything other than visualization? If you change the visualization projection (shifting the meridian (+lon_0) for worldwide distributions) repeatedly, you will have to clip to suit. If visualization, does it make sense to modify the meridian (the reader may see the shifts as confusing)? You are trying to deal with multiple somethings that you are interpolating (to a grid?) and displaying in possibly different windows possibly depending on the extent of the prediction grid. Your workflow isn't obvious. You haven't explained why you need the projection so far (other than aesthetics, maybe?). You may be able to interpolate on the sphere, so unless the target grid cells must be of fixed area - and possibly not even then, you don't need to make life complicated. There are problems in creating global visualizations in any case, as the lines/polygons need clipping at +/- 180 degrees from the meridian. Roger Dr. Francesco Carotenuto (Ph.D.) Department of Earth Sciences Federico II University Largo S. Marcellino 10, 80138 NAPLES (ITALY) http://francesco-carotenuto.blogspot.it/ http://3dpaleontology.blogspot.it/ Il Giovedì 5 Dicembre 2013 21:23, Roger Bivand roger.biv...@nhh.no ha scritto: On Thu, 5 Dec 2013, Francesco Carotenuto wrote: Dear list, as in the title, I found that when projecting a world map by mollweide projection and setting a?? longitude value, some strange parallels are drawn. The real problem is that these lines ???interact??? with a raster when performing the ???over??? function. How can I manage it? The underlying question is what you want to do with the objects after projection. The problem in your example is that some polygons/lines cross -180+66 degrees, so get reflected back across. A possible solution is to clip the polygons before projection at -180+66 degrees +/- a small number, but before looking at this in more detail, I need to know what the projected object is intended for. When +lon_0=0, there is no such effect, as this data set is clipped at -180 already. Roger Here I provide a simulated example to let you understand what I mean Thanks.?? #Example of a world map mollweide projection with a specific longitude value library(sp) library(maptools) library(rgdal) data(wrld_simpl) w_pol-wrld_simpl w_pol-spTransform(w_pol, CRS(+proj=moll +lon_0=66)) plot(w_pol) ?? # Intersecting a grid with the projected world map gr-expand.grid(x = seq(-180, 180, length.out = 100), y = seq(-90, 90, length.out = 100)) coordinates(gr)-~x+y proj4string(gr)-CRS(+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs) gr-spTransform(gr, CRS(+proj=moll +lon_0=66)) plot(gr[w_pol,]) ?? ?? Dr. Francesco Carotenuto (Ph.D.) Department of Earth Sciences Federico II University Largo S. Marcellino 10, 80138 NAPLES (ITALY) http://francesco-carotenuto.blogspot.it/ http://3dpaleontology.blogspot.it/?? [[alternative HTML version deleted]] -- Roger Bivand Department of Economics, Norwegian School of Economics, Helleveien 30, N-5045 Bergen, Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43 e-mail: roger.biv...@nhh.no -- Roger Bivand Department of Economics, Norwegian School of Economics, Helleveien 30, N-5045 Bergen, Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43 e-mail: roger.biv...@nhh.no ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] 3rd+ rasters added using plot(r, add=TRUE) are misregistered
Hi Michael, Thanks for those thoughts. The issue you raise (which has bit me before as well) appears to be unrelated. In particular, your resizing issue disappears if one sets useRaster=FALSE, whereas my problem persists whether or not rasterImage() is used to draw the raster. Try this to get a version of the volcano and contour that resize/stretch in sync: plot(raster(v), xaxs = i, useRaster=FALSE) lines(cl) ## Check that resizing works correctly. (It does, presumably b/c ## image(..., useRaster=FALSE) doesn't attempt to maintain aspect ratio of the raster. The problem I described seems related to resetting of elements of `par` by the first call to `plot(r, ..., add=TRUE)`, as can be seen by doing: plot(r, col=yellow, legend=FALSE) p0 - par() plot(r, col=green, legend=FALSE, add=TRUE) p1 - par() plot(r, col=grey, legend=FALSE, add=TRUE) p2 - par() plot(r, col=yellow, legend=FALSE, add=TRUE) p3 - par() ## Find elements of par that are changed by first call to plot(r, ..., add=TRUE) names(p1)[!mapply(identical, p0,p1)] # [1] cxy pin plt ## Show that par does not change thereafter all(identical(p1,p2), identical(p2,p3)) # [1] TRUE As might be expected if par is to blame, frequently resetting par to its 'original' value makes the plots come out right: plot(r, col=yellow, legend=FALSE) p0 - par(no.readonly=TRUE) plot(r, col=green, legend=FALSE, add=TRUE) plot(r, col=blue, legend=FALSE, add=TRUE) par(p0) plot(r, col=red, legend=FALSE, add=TRUE) plot(r, col=grey, legend=FALSE, add=TRUE) par(p0) plot(r, col=dodgerblue, legend=FALSE, add=TRUE) The piece that is most puzzling to me is that, although cxy, pin, plt get reset after the first added layer is plotted, it isn't until _two_ more layers are added that the problem manifests itself Seems like an interesting puzzle here, which might be solvable by doing `debug(raster:::.imageplot)` and stepping through it, paying careful attention to how par gets set and reset during each of the calls to plot(). Cheers, Josh On Thu, Dec 5, 2013 at 2:02 PM, Michael Sumner mdsum...@gmail.com wrote: I just wonder if this is related to an existing resize problem when overplotting? I see this on the windows() graphics device, you need to try these two sets of plot commands and resize the window (especially stretch it in one direction) to see the problem at ## 1 and ## 2 v - list(x = seq(0, 1, length = nrow(volcano)), y = seq(0, 1, length = ncol(volcano)), z = volcano) cl - contourLines(v, levels = 160)[[2]] image(v, xaxs = i) lines(cl) ## (1) now resize the window, it's all good - the contour stays in the ## right place ## however, resize after adding the lines to a raster plot ## and the relationship is broken library(raster) plot(raster(v), xaxs = i) lines(cl) ## (2) May be totally unrelated but I thought it worth pointing out. It's hard to trace through the setting up that raster's plot does (well, I find it hard to do - I've tried a few times to figure out where these problems are happening). Cheers, Mike. On Wed, Dec 4, 2013 at 5:57 AM, Josh O'Brien joshmobr...@gmail.com wrote: Hi all, Briefly, running the following code will not produce a single grey ring, as would be expected. The first three raster layers all stack up neatly, but the 4th and subsequent raster layers are all shifted/skewed to the right. It appears that the plotting engine loses track of the fact that some space in the original plot was allocated to a legend, and just begins plotting as if a wider canvas were available to it. I'm not sure whether this is ultimately an issue with raster or sp or something else, which is why I'm noting it here on R-sig-geo. (I'm also fully aware that there are other and better ways to add several raster layers to a plot, but still think the following behavior qualifies as a bug that I ought to report.) library(maptools) ## Only needs to be installed for example data library(raster) library(rgeos) ## Create an example raster p - shapefile(system.file(shapes/co37_d90.shp, package=maptools)) p - p[31,] ## A tall narrow county polygon pr - gDifference(gBuffer(p, width=.01), p) r - rasterize(pr, raster(extent(pr), ncol=100, nrow=100)) ## These three are properly registered on one another plot(r, col=yellow, legend=FALSE) plot(r, col=green, legend=FALSE, add=TRUE) plot(r, col=grey, legend=FALSE, add=TRUE) ## All subsequent layers are improperly shifted/skewed to right plot(r, col=yellow, legend=FALSE, add=TRUE) plot(r, col=blue, legend=FALSE, add=TRUE) plot(r, col=red, legend=FALSE, add=TRUE) plot(r, col=grey20, legend=FALSE, add=TRUE) ## Following the above, SpatialPolygons are also shifted/skewed plot(p) My session info, in case it matters, is: sessionInfo() R version 3.0.2 (2013-09-25)
Re: [R-sig-Geo] 3rd+ rasters added using plot(r, add=TRUE) are misregistered
Interesting and difficult! Thanks for pointing out the useRaster difference, I had mixed up TRUE/FALSE in my original report (offlist) on that one. That might help me find the solution, though sorry I have nothing for your issue. :) Cheers, Mike. On Fri, Dec 6, 2013 at 11:58 AM, Josh O'Brien joshmobr...@gmail.com wrote: Hi Michael, Thanks for those thoughts. The issue you raise (which has bit me before as well) appears to be unrelated. In particular, your resizing issue disappears if one sets useRaster=FALSE, whereas my problem persists whether or not rasterImage() is used to draw the raster. Try this to get a version of the volcano and contour that resize/stretch in sync: plot(raster(v), xaxs = i, useRaster=FALSE) lines(cl) ## Check that resizing works correctly. (It does, presumably b/c ## image(..., useRaster=FALSE) doesn't attempt to maintain aspect ratio of the raster. The problem I described seems related to resetting of elements of `par` by the first call to `plot(r, ..., add=TRUE)`, as can be seen by doing: plot(r, col=yellow, legend=FALSE) p0 - par() plot(r, col=green, legend=FALSE, add=TRUE) p1 - par() plot(r, col=grey, legend=FALSE, add=TRUE) p2 - par() plot(r, col=yellow, legend=FALSE, add=TRUE) p3 - par() ## Find elements of par that are changed by first call to plot(r, ..., add=TRUE) names(p1)[!mapply(identical, p0,p1)] # [1] cxy pin plt ## Show that par does not change thereafter all(identical(p1,p2), identical(p2,p3)) # [1] TRUE As might be expected if par is to blame, frequently resetting par to its 'original' value makes the plots come out right: plot(r, col=yellow, legend=FALSE) p0 - par(no.readonly=TRUE) plot(r, col=green, legend=FALSE, add=TRUE) plot(r, col=blue, legend=FALSE, add=TRUE) par(p0) plot(r, col=red, legend=FALSE, add=TRUE) plot(r, col=grey, legend=FALSE, add=TRUE) par(p0) plot(r, col=dodgerblue, legend=FALSE, add=TRUE) The piece that is most puzzling to me is that, although cxy, pin, plt get reset after the first added layer is plotted, it isn't until _two_ more layers are added that the problem manifests itself Seems like an interesting puzzle here, which might be solvable by doing `debug(raster:::.imageplot)` and stepping through it, paying careful attention to how par gets set and reset during each of the calls to plot(). Cheers, Josh On Thu, Dec 5, 2013 at 2:02 PM, Michael Sumner mdsum...@gmail.com wrote: I just wonder if this is related to an existing resize problem when overplotting? I see this on the windows() graphics device, you need to try these two sets of plot commands and resize the window (especially stretch it in one direction) to see the problem at ## 1 and ## 2 v - list(x = seq(0, 1, length = nrow(volcano)), y = seq(0, 1, length = ncol(volcano)), z = volcano) cl - contourLines(v, levels = 160)[[2]] image(v, xaxs = i) lines(cl) ## (1) now resize the window, it's all good - the contour stays in the ## right place ## however, resize after adding the lines to a raster plot ## and the relationship is broken library(raster) plot(raster(v), xaxs = i) lines(cl) ## (2) May be totally unrelated but I thought it worth pointing out. It's hard to trace through the setting up that raster's plot does (well, I find it hard to do - I've tried a few times to figure out where these problems are happening). Cheers, Mike. On Wed, Dec 4, 2013 at 5:57 AM, Josh O'Brien joshmobr...@gmail.com wrote: Hi all, Briefly, running the following code will not produce a single grey ring, as would be expected. The first three raster layers all stack up neatly, but the 4th and subsequent raster layers are all shifted/skewed to the right. It appears that the plotting engine loses track of the fact that some space in the original plot was allocated to a legend, and just begins plotting as if a wider canvas were available to it. I'm not sure whether this is ultimately an issue with raster or sp or something else, which is why I'm noting it here on R-sig-geo. (I'm also fully aware that there are other and better ways to add several raster layers to a plot, but still think the following behavior qualifies as a bug that I ought to report.) library(maptools) ## Only needs to be installed for example data library(raster) library(rgeos) ## Create an example raster p - shapefile(system.file(shapes/co37_d90.shp, package=maptools)) p - p[31,] ## A tall narrow county polygon pr - gDifference(gBuffer(p, width=.01), p) r - rasterize(pr, raster(extent(pr), ncol=100, nrow=100)) ## These three are properly registered on one another plot(r, col=yellow, legend=FALSE) plot(r, col=green, legend=FALSE, add=TRUE) plot(r, col=grey, legend=FALSE, add=TRUE) ## All subsequent layers are improperly
[R-sig-Geo] original coordinates of the redwood data?
Lee De Cola ldec...@comcast.net writes: does anyone know the original, true coordinates of the spatstat redwoodfull data? i can't find them in: Strauss, D. J. (1975). A Model for Clustering. Biometrika 62(2): 467-475. i think reporting spatial data in rescaled units is unscientific. Please first read the help file for the dataset. (I think it is unprofessional to send out a question before you have read the documentation %^]) The help file says that these data were scanned manually from a reprint of Strauss's paper and that the dimensions of the figure are only *approximately* known to be about 128 feet across. The paper by Strauss does not give the exact dimensions of the figure, and does not give the original coordinates. The data are lost, according to David Strauss. Hence we were obliged to scan a reprint - which has the risk of inaccuracy (coincident points may be lost; the x and y axes may have been unequally scaled, nonlinearly scaled, etc). If you are confident enough about the dimensions, it is easy to rescale the data: X - rescale(redwoodfull, 128) unitname(X) - c(foot, feet) In spatstat we do not attribute physical scales to datasets unless the scale is reliable, and that is why 'redwoodfull' is not scaled to 128 feet. i think reporting spatial data in rescaled units is unscientific. Sure, the redwood dataset does not meet modern standards of data integrity. The redwood patterm is a 'legacy' dataset that has been widely used and re-used over the years to compare different methods (a bit like Fisher's iris data or the sunspot series). That is why it is included in spatstat. In the 1970's it was actually standard practice to rescale point patterns to the unit square, to avoid numerical difficulties (for example all the data in Peter Diggle's 1983 book are rescaled to the unit square). Also it was a lot harder to obtain data, and scientists were even more protective of their data than they are now, so spatial data were not widely shared, and when shared it was often 'anonymised' by removing the spatial scale so that competitors could not re-analyse the data. Adrian Baddeley author of 'spatstat' Prof Adrian Baddeley FAA University of Western Australia ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
