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

[Roger Bivand]

On Thu, 11 Jan 2007, David Forrest wrote:

> Hi All,
> 
> I'd like to make graphic files of plots with controlled pixel->coordinate 
> mapping.  Ultimately, I'd like to get the graphs out into Google earth 
> KML, but to to that, I need good control of the margins or framing around 
> the plot.
> 
> Basically, I'm using a spatial polygons dataframe with some 
> finite element (30K elements) model output.
> 
> png(file='inun.png',width=2000,height=2000,bg = "transparent")
> plot(cb,xaxs='i',yaxs='i',
>     xlim=bbox(cb)[1,],ylim=bbox(cb)[2,],lty=0,col=cb$inun*4)
> dev.off();
> 
> but I still have an uneven margin/border around my graphic that is hard to 
> correlate to the bbox.

Could you try rasterising your output to a raster in geographical
coordinates, creating the RGB bands by hand, and using functions in the
rgdal package to emit the PNG and WLD files? This isn't your case, but has
worked for the Meuse Bank data set.

The key bits are defining a bounding mask (if your data are not 
rectangular in geographical coordinates), creating a GridTopology in 
geographical coordinates, making a SpatialGrid, adding the mask and 
reducing to SpatialPixels, for you overlaying the SpatialPolygonDataFrame 
on the SpatialPixels to get a data class value for each pixel, put the 
data class values in the SpatialPixelsDataFrame, use vec2RGB with your 
choice of break points and palette to generate an RGB matrix, add red, 
green, and blue columns to the SpatialPixelsDataFrame:

RGB <- vec2RGB(llSPix$pred, breaks=fj7$brks, col=rev(bpy.colors(7)))
llSPix$red <- RGB[,1]
llSPix$green <- RGB[,2]
llSPix$blue <- RGB[,3]
fullgrid(llSPix) <- TRUE
llSPix$alpha <- as.integer(ifelse(is.na(llSPix$pred), 0, 255))

and finally use the alpha channel to drop the out-of-mask pixels after 
going back to a SpatialGridDataFrame. Having to set the three colour 
bands manually is rather clunky, but isn't that what scripts are for? Then 
write out:

dx <- create2GDAL(llSPix[c("red", "green", "blue", "alpha")], 
drivername="GTiff", type="Byte", options="ALPHA=YES")
png_dx <- copyDataset(dx, driver="PNG", options="WORLDFILE=YES")
GDAL.close(dx)
saveDataset(png_dx, filename="log_zinc.png", options="WORLDFILE=YES")
GDAL.close(png_dx)

The round-the-houses is needed because the GDAL PNG driver cannot create 
directly, just by copying, so we write to geotiff with alpha support and 
copy across to PNG with a world file. Open in GE adjusting the opacity to 
suit your taste, probably entering the bounds into the location tab by 
hand. The values from gdalinfo of the *.png file are fine, or the bbox() 
of the SpatialGrid object. If anyone knows how to get GE to use a *.wld 
file and read the location correctly, please say!

We can iterate this if you like, but once you've rasterised, the rest is 
feasible. With the output raster in geographical coordinates, it seems to 
work OK.

Best wishes,

Roger

> 
> Any hints?
> 
> 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