Skip to content

spplot with two rasters

14 messages · Edzer Pebesma, Thomas Adams, Dylan Beaudette +3 more

Original
Dylan Beaudette
# 
Hi,

Is it possible to plot two raster images using spplot() in a manner similar 
to:

pts <- list("sp.points", points_file, pch = 4, col = "black", cex=0.5)
spplot(raster_file, zcol="elev.pred", sp.layout=list(pts))

Note that one of the raster images is an aerial photo, used only for context, 
while the second one is one with interesting z-values. The second raster is 
masked and thus does not cover the entire region.

Thanks,

Dylan
Edzer Pebesma
# 
I find it hard to imagine how you want to plot two raster files on top 
of each other. Do you want some form of transparency? If it is just one 
overlaying the other, you could use overlay to find out which cells in 
raster 1 to replace with those in raster 2 before plotting.
--
Edzer
Dylan Beaudette wrote:
Dylan Beaudette
#
On Tuesday 04 March 2008, Edzer Pebesma wrote:
Hi Edzer,

I generally agree that plotting one raster file "over" another raster file 
would be of little use. In this case, one of the raster files (the 
interesting one) has been masked with nodata, such that it only really covers 
about 30% of the region of interest. The other raster is just contextual 
data, and thus would be useful to plot "behind" the first raster.

Ideas?

Dylan

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Thomas Adams
    

    

      
      #
    

  


  

      
Dylan,

I think a solution using GRASS can be found on pages 110-11 of "Open 
Source GIS: A GRASS GIS Approach", 3rd Ed. The same material is covered 
in the 2nd Ed. as well, where you use r.mapcalc to combine two rasters 
and judicious use of MASKs; a conditional statement in r.mapcalc is the key.

Regards,
Tom

        
Dylan Beaudette wrote:

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Thomas Adams
    

    

      
      #
    

  


  

      
Dylan,

I think a solution using GRASS can be found on pages 110-111 of "Open 
Source GIS: A GRASS GIS Approach", 3rd Ed. The same material is covered 
in the 2nd Ed. as well, where you use r.mapcalc to combine two rasters 
and judicious use of MASKs; a conditional statement in r.mapcalc is the key.

Regards,
Tom

        
Dylan Beaudette wrote:

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    
6 days later

    

      
      

        


  

        

      Dylan Beaudette
    

    

      
      #
    

  


  

      
On Thursday 06 March 2008, Thomas Adams wrote:

            

      
      
      
    

        
Hi Tom,

Thanks for the suggestion. This works, but we were hoping to make the plot in 
R-- as the PDF output is hard to match with GRASS alone. I suppose I will 
just have to try using two rasters with spplot() and see what happens. 

Cheers,

Dylan

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Roger Bivand
    

    

      
      #
    

  


  

      
On Wed, 12 Mar 2008, Dylan Beaudette wrote:

        

            

      
      
      
    

        
The problem is that customising a panel to insert an image probably means 
writing just that, customised panels, and spplot() essentially does that 
already. You could look at Chapter 5 in Paul Murrell's book, to see how to 
insert grid output into lattice output (to try to put an image as a 
backcloth, but there is no example for this case - the example is for 
adding a location position.

My guess would be that if you are only displaying a single panel, you will 
find it easier to use base graphics, and simply say image() twice, and use 
legend() (or legend.krige() from geoR).

Roger

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Agustin Lobo
    

    

      
      #
    

  


  

      
Well, the problem is that R does not have a real geographical
display. While things can be done going back and forth from R to GIS,
this procedure soon becomes very inconvenient. It's ok
for learning and teaching, but not for real applications.
Maybe getting an existing GIS to display spatial R objects
is actually easier than developing a geographical display for R.
Agus

Dylan Beaudette escribi?:

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Edzer Pebesma
    

    

      
      #
    

  


  

      
Agus, I disagree with your claim that R plot functions are not ok for 
real applications. Also, I use real applications for learning and 
teaching -- they are not fundamentally different.

In the trellis panel functions, or in the direct plot/image/lines etc 
functions, you plot in data coordinates. The main thing that the plot 
and spplot methods in package sp control is the aspect ratio. How else 
do data coordinates differ from geographical coordinates when it comes 
to plotting?

Are you refering to taking care of geographic (re)projection at that level?
--
Edzer

        
Agustin Lobo wrote:

            

      
      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Roger Bivand
    

    

      
      #
    

  


  

      
On Thu, 13 Mar 2008, Edzer Pebesma wrote:

        

            

      
      
      
    

        
Exactly. In fact I think that a grid panel might make it possible to 
insert an image backdrop behind spplot(), which was the original question 
from Dylan - it would make a nice student project in visualisation. 
Something like an sp.image() for sp.layout=, but using its own palette and 
being painted over by the "real" data.

Roger

            

      
      
      
    

        

      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Agustin Lobo
    

    

      
      #
    

  


  

      
Edzer,

Edzer Pebesma escribi?:

            

      
      
      
    

        
Yes they are: we use the results of real applications as examples
for teaching. Real applications need a lot of (geo)graphical interaction
for being developed, which cannot be done in the (not geo)graphical
display of R, you need a GIS display for that. In my case, in most cases 
I have to export the
sp objects as rasters or shp files, bring them to a GIS and interact
with them (and many other layers of information) there (for example,
getting on-screen values as tips, zooming in and out, changing the
background from one color composite to another, calculating buffers,
interactively measuring distances, selecting polygons according to an
associated database...). These operations are required to evaluate
results and to think on alternative approaches, and have to be done
easily and fast. If your problem is
solved with one single cycle involving one analysis in R, exporting to shp
and displaying, this is fine. But if you have to iterate
this process, having the option of displaying the sp
objects directly in the GIS display would be an advantage.

            

      
      
      
    

        
This is one case, yes: if your GIS supports on-the-fly reprojection, you
can use layers from geodata distribution centers that work with a
different projection
than yours and display them along with your R (exported) layers,
no need for the user to run a reprojection, it's automatically done.
This is often the case here, where we normally work with UTM but get
many layers at European scale that are in a different projection.

In a previous message you said:

            

      
      
      
    

        
Well, in many GIS displays you just define a rectangle that you can move
on top of one raster as a window through which you can see (and check
values of) an underlying raster. This is most useful, and is set up in
seconds. Also, perhaps Dylan's problem could have been solved in a
GIS by
displaying the Z values as terrain, and draping the other raster on top
of it. Or just defining shades for the second raster. Or perhaps setting
a 50% of transparency for some regions of one raster...

Geographic information requires GIS display, that's why GIS displays exist.

In much the same way as the analysis of geographic data required the
design of sp
objects, the visualization of geographic data requires a geographic
display. Getting sp objects directly displayed on a GIS would be great.

Agus

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Hadley Wickham
    

    

      
      #
    

  


  

          

      
      
      
    

        
Another approach would be to develop a better interface between
spatial data and ggplot2, which provides much better support for layer
compositing (and in general for modifying a plot after you have
created it).  There are some simple examples at
http://had.co.nz/ggplot2/geom_polygon.html and
http://had.co.nz/ggplot2/coord_map.html, but I don't know enough about
the problems of spatial data to be able to do it all myself.

Hadley

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Roger Bivand
    

    

      
      #
    

  


  

      
On Thu, 13 Mar 2008, hadley wickham wrote:

        

            

      
      
      
    

        
Maybe, but as Wilkinson says, his proposal is not a GIS either. I'm sorry 
nobody replied to your question earlier about using mapproj. Even the way 
Wilkinson addresses things isn't going to deliver unless the data and the 
graphics stay very close together.

This is the mapproj problem, the projected coordinates are an eclectic 
rendering, not something that can be used for anything else than display. 
Over thousands of km, the problems are not so clear, but over a couple of 
km they are showstoppers.

It looks as though Simon Urbanek has been trying to "do it his way" for 
iplots by publishing a limited interface to PROJ.4 on CRAN as a proj4 
package, presumably because there are real issues in providing proper 
support on OSX, and on Windows (I don't imagine that he didn't know that 
PROJ.4 is interfaced in rgdal, because I have asked for his help in 
finding a mechanism for producing rgdal binaries for OSX with minimal 
intervention). The Windows issues are fully addressed by rgdal (but not in 
proj4, there is no PROJ.4 share/), which he could have used, but chose not 
to. Perhaps I should ask again whether he could help automate the 
production of rgdal binaries for OSX?

Yes, you are right that this needs doing in concert between people who 
understand both sides, the graphics side, and crucially the spatial data 
side, because the layering fails unless the data are correctly aligned. 
This constrains the ways the data can be represented, and the S data model 
of a two-column matrix with NAs to separate geometric entities is simply 
not sufficiently robust.

Roger

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    
3 days later

    

      
      

        


  

        

      Dylan Beaudette
    

    

      
      #
    

  


  

      
On Wednesday 12 March 2008, Roger Bivand wrote:

            

      
      
      
    

        
Thanks for the tips Roger. I think that in this case it will be simplest to 
use multiple calls to image().

Cheers,

Dylan