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

[David Forrest]

On Wed, 17 Jan 2007, Roger Bivand wrote:

> 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.


I've solved this somewhat by explicitly setting the xlim and ylim.  With 
Using the theme.novpadding from above, this seems to work:
library(maptools)
library(sp)
library(stats)

qk<-SpatialPointsDataFrame(quakes[,c(2:1)],quakes)
proj4string(qk)<-CRS("+proj=longlat")
qk.p<-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),
        xlim=bbox(qk)[1,],ylim=bbox(qk)[2,])

print(qk)

# With a controlled-size png file:
myheight=1024;mywidth=myheight/mapasp(qk)
png(file='mymap.png',width=mywidth,height=myheight)
print(qk.p)
dev.off()
bbox(qk)

However, now I'd like to be able to eliminate the key entirely.  It seems 
that the keys are handled differently between points, grid, lines, and 
polygons, and I'm not sure there is a general control at the spplot() 
level.

> 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

Exactly.

> 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.

Yes.  I do think that may be viable as well, but it does depend on making 
sure the plot does not draw any extraneous info like colorbars or keys. 
I'm not terribly skilled with OO programming and digging out just what 
happens behind a spplot(SpatialObject,...) or plot(SpatialObject) seems 
confusing.  (Besides 'plot' 'spplot', or '?plot' or '?ssplot' how do you 
see what happens next?)

Alternately, it might be enough to somehow report the extents of the 
graphics device in data coordinates.

> 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.

Hmm.  I'll look at that.

> 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
>>
>
>

-- 
  Dr. David Forrest
  drf at vims.edu                                    (804)684-7900w
  drf5n at maplepark.com                             (804)642-0662h
                                    http://maplepark.com/~drf5n/