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

[David Forrest]

On Thu, 11 Jan 2007, Roger Bivand wrote:

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

At the core, my finite element model output is defined on either the 
elements (triangular polygons on a unstructured mesh) or nodes (points) 
and I'm looking for a way to create the graphics in different usable ways. 
I have a shapefile of my finite element grid in geographical coordinates, 
but I have a long multivariate timeseries for the element (polygons) that 
doesn't translate well to shapefiles.  I use R (and sometimes Matlab) for 
visualization.

Doing this gets me good graphics for day-to-day visualization:

cb<-readShapePoly('SHAP/grid.shp')
xx<-read.table('output.txt')
cb$inun<-xx[cb$ID,2]  #1# an indicator variable for inundation
plot(cp)              #  shows the 30k Triangle grid 
plot(cb,col=cb$inun*4,lty=0)  # Shows the inundated cells

Since the data is a multivariate timeseries, I've been saving captures of 
the graphics for animations of time and other explorations of the model 
data by mostly varying from the step #1# above, with many rasterizings.

For right now, it looks like I've got an old rdgal without vec2RGB, and 
might have to manually georeference plots from my static domain.

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

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