Skip to content

adding points to spplot grid

5 messages · Edzer Pebesma, Paul Hiemstra, Ben Madin +1 more

Original
Nick Matzke
# 
Hi,

I have been working with spplot.  I can get continuous grids 
to plot with spplot, but what do I do if I want to plot some 
points with x,y coordinates on top of the image?

Cheers!
Nick
Edzer Pebesma
#
On 11/06/2010 01:37 AM, Nick Matzke wrote:
http://r-spatial.sourceforge.net/gallery/
Paul Hiemstra
#
On 11/06/2010 01:37 AM, Nick Matzke wrote:
?spplot, look at the sp.layout argument.

cheers,
Paul
Ben Madin
#
On 06/11/2010, at 9:03 PM, Edzer Pebesma wrote:

            
        
Nice resource, thanks for the heads up.

cheers

Ben

            

      
      
      
    

        

      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Nick Matzke
    

    

      
      #
    

  


  

      
On 11/6/10 3:03 AM, Edzer Pebesma wrote:

        
>
 >

        
> On 11/06/2010 01:37 AM, Nick Matzke wrote:

        
>> Hi,
 >>
 >> I have been working with spplot.  I can get continuous 
grids to plot
 >> with spplot, but what do I do if I want to plot some 
points with x,y
 >> coordinates on top of the image?
 >>
 >> Cheers!
 >> Nick
 >>
 >>
 > http://r-spatial.sourceforge.net/gallery/


Thanks! This helped although it took some digging.  For 
others with this question:

=================================
# I start with some (randomly generated for now)
# * TRUE/FALSE values (samp_hits)
# * x & y columns giving sample locations (sample_locs)

# generate the tmpdata for species #1 (samp_hits, column #1)
tmpdata = as.data.frame(cbind(sample_locs[,1], 
sample_locs[,2], samp_hits[,1]))

names(tmpdata) = c("x", "y", "sp1")
coordinates(tmpdata) = ~x+y
# set projection to lon (x) & lat (y), i.e. unprojected;
# find projections with Terminal: proj -lp
proj4string(tmpdata) = CRS("lonlat")

hits = tmpdata[tmpdata$sp1 == 1, ]
misses = tmpdata[tmpdata$sp1 == 0, ]


#plot points for kicks
plot(0,0, xlim=minmax_literal(sample_locs[, 1]), 
ylim=minmax_literal(sample_locs[, 2]) )

points(hits, col="blue", pch=19)
points(misses, col="blue", pch=1)


# Make variogram

default_vgm = vgm(0.15, "Exp", 80, 0)
gs = gstat(id="sp001", formula=sp1~1, data=tmpdata, beta = 
0.1,  model=default_vgm)
vg = variogram(gs)
plot(vg)


# Make a tmpgrid
xs=xmin:xmax
ys=ymin:ymax

# use expand.grid to make all possible combinations of inputs
tmpgrid = expand.grid(xs, ys)
names(tmpgrid) = c("x", "y")
gridded(tmpgrid) = ~x+y
# set projection to lon (x) & lat (y), i.e. unprojected;
# find projections with Terminal: proj -lp
proj4string(tmpgrid) = CRS("lonlat")



# Fit a variogram model to the empirical variogram
# for single variogram
sp1_vgm = fit.variogram(vg, model=default_vgm)

plot(vg, model=default_vgm)
plot(vg, model=sp1_vgm)

# simple kriging
sk = krige(sp1~1, tmpdata, tmpgrid, sp1_vgm, nsim=0, 
indicators = TRUE, nmax = 40)

# plot the result
print(spplot(sk, cuts=1, col.regions=grey.colors(2, 0.6, 
0.90, 2.2)))
print(spplot(sk))


alphaChannelSupported = function()
	{
	!is.na(match(names(dev.cur()), c("pdf")))
	}

spplot_hits = list("sp.points", hits, pch=19, col="black", 
alpha = ifelse(alphaChannelSupported(), .5, 1))
spplot_misses = list("sp.points", misses, pch=21, cex=0.75, 
col="black", alpha = ifelse(alphaChannelSupported(), .5, 1))

# col.regions does color scale for grid
# sp.layout contains a list of the various points, text, 
etc. to plot
# at=seq() does limits of the color bar
spplot(sk["var1.pred"], scales=list(draw = TRUE), 
col.regions=rev(rainbow(100)[1:75]), 
sp.layout=list(spplot_hits, spplot_misses), at=seq(0, 1, 1/75))

=================================