[R-sig-Geo] Marginless plot output for georegistration of output graphics

[Roger Bivand]

On Tue, 16 Jan 2007, David Forrest wrote:

> On Thu, 11 Jan 2007, Roger Bivand wrote:
> 
> > On Thu, 11 Jan 2007, David Forrest wrote:
> >
> >> Hi All,
> ...
> 
> I think I can make marginless plots in base, lattice, and sp graphics now, 
> but I'm still missing something about the data ranges for spplot.
> 
> #################################
> # Marginless plots in base graphics:
> library(stats)
> par.save<-par(no.readonly = TRUE);
> par(mar=c(0,0,0,0))
> plot(quakes$long,quakes$lat,xaxs='i',yaxs='i',xlab=FALSE,axes=FALSE)
> par(par.save)
> 
> ########### in lattice/trellis/grid
> library(lattice)
> theme.novpadding <-
> # From Deepayan Sarkar's posting of 2006-10-12:
>     list(layout.heights =
>          list(top.padding = 0,
>              main.key.padding = 0,
>              key.axis.padding = 0,
>              axis.xlab.padding = 0,
>              xlab.key.padding = 0,
>              key.sub.padding = 0,
>              bottom.padding = 0),
>          layout.widths =
>          list(left.padding = 0,
>              key.ylab.padding = 0,
>              ylab.axis.padding = 0,
>              axis.key.padding = 0,
>              right.padding = 0))
> print( # lattice/grid With device extents == data extents
> xyplot(lat ~ long , data=quakes, pch=".",
>         ,xlab=NULL,ylab=NULL
>         ,scales=list(axs='i',draw=FALSE)
>         ,par.settings=theme.novpadding))
> 
> #########  In spatial (with internal & device margin)
> # An example of a marginless plot in library(sp)
> library(sp)
> qk<-SpatialPointsDataFrame(quakes[,c(2:1)],quakes)
> proj4string(qk)<-CRS("+proj=longlat")
> spplot(qk,'mag',key.space=list(x=0.25,y=.2,corner=c(1,1)),
>         par.settings=theme.novpadding,scales=list(axs='i',draw=FALSE))
> 
> #####################
> 
> The device/aspect ratio margin is controllable and manageable with a 
> mapasp(gk) based device call.
> 
> I'd still like to eliminate the remaining internal margin from around the 
> plot.  I saw it once in my experimenting, but I am not sure if it was some 
> interaction with par() or something and have not able to duplicate 
> it.

What you are trying to do is effectively to implement a graphics device 
that only takes the printed area. For many polygons, I've tried overlaying 
on a SpatialGrid, but this is too time-consuming. I also looked at 
plotting to postscript, resetting the bounding box to be as tight as 
possible (script bbox.pl originally by Dov Grobgeld using GhostScript), 
and converting (Image Magic or other) to PNG, but am not there yet, 
although I think it is viable. The eps output is in its own coordinates, 
though, so the bounding box for GE will need to be extracted first. 

Another possibility is the xfig format and fig2dev, and looks promising:

library(maptools)
xx <- readShapePoly(system.file("shapes/sids.shp", package="maptools")[1], 
  IDvar="FIPSNO", proj4string=CRS("+proj=longlat +ellps=clrk66"))
tf <- tempfile()
xfig(tf)
plot(xx)
dev.off()
tf2 <- tempfile()
system(paste("fig2dev -L png", tf, tf2))
system(paste("display", tf2))

does seem to work - using the width and height arguments in xfig() can add 
pixels, and extra arguments to fig2dev can change things too.

If anyone has lots of time to spare, a SpatialGridDataFrame device with 
red, green, blue and alpha channels would be nice (or a GDAL MEM device, 
or whatever).

Roger

> 
> Dave
> 

-- 
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no