[R-sig-Geo] Vector to raster?

[Roger Bivand]

On Sat, 14 Feb 2009, Jonathan Greenberg wrote:

> Thanks -- another few questions along these lines:

Perhaps you could simplify by providing code examples with a built-in data 
set and so on? Verbatim copies do let others see what is going on. So I'll 
start by trying to reconstruct what you seem to be asking about:

library(rgdal)
dsn <- system.file("vectors", package = "rgdal")[1]
layer <- "scot_BNG"
scot_BNG <- readOGR(dsn, layer)

library(maptools)
SG <- Sobj_SpatialGrid(scot_BNG)$SG
o <- overlay(scot_BNG, SG)
SGDF <- SpatialGridDataFrame(slot(SG, "grid"),
  proj4string=CRS(proj4string(SG)), data=o)
spplot(SGDF, "SMR", col.regions=bpy.colors(20))

The Sobj_SpatialGrid() function in maptools takes a maxDim= argument, 
which indirectly controls the (square) cell resolution. It is used in 
creating PNG devices for displaying arbitrary objects in geographical 
coordinates on Google Earth, so is related to vector to raster. An 
alternative is to create a GridTopology object from scratch, here with 
10km by 10km cells:

bb <- bbox(scot_BNG)
grd <- GridTopology(cellcentre.offset=floor(bb[,1]), cellsize=c(10000,
  10000), cells.dim=ceiling(diff(t(bb))/10000))
o <- overlay(scot_BNG, SpatialGrid(grd))
SGDF <- SpatialGridDataFrame(grd,  proj4string=CRS(proj4string(scot_BNG)),
  data=o)
spplot(SGDF, "SMR", col.regions=bpy.colors(20))

>
> 1) Is there any way to determine the type of vector file in advance of 
> defining a layer name?  How can I list layer names affiliated with a given 
> vector (it seems to be part of ogrinfo from the main gdal release, but when i 
> do a ogrinfo(filename) I get yelled at to feed it a layer name first).

ogrInfo(dsn, layer)
system(paste("ogrinfo", dsn))
# only works if ogrinfo is in your PATH
list.files(dsn)

suggests that your description isn't adequate - that is my external 
ogrinfo program appears only to choose one driver, based on the fact that 
it has been passed a directory. If the dsn was suited to other drivers, it 
would have found them, but so would, for example:

list.files(dsn, pattern="shp$")

> 2) Is there a way to tell the feature type (point, polygon, line, etc...) of 
> a readOGR object?

l <- list.files(dsn, pattern="shp$")
sapply(l, function(fn) getinfo.shape(paste(dsn, fn, sep="/")))

for shapefiles and using a function in maptools, type is as in the ESRI 
documentation. I'll try to add a similar facility to rgdal for available 
drivers - in fact, the type is declared for each feature, but is assumed 
uniform on return after accessing all features.

> 3) How do I get a polygon/line's node coordinates?  Using 
> coordinate(readOGR(...)) I only get a single coordinates associated with each 
> polygon/line in the vector.

What do you mean by "node coordinates" - no topology is being done here, 
so no arc/node analysis is available. If you mean the coordinates from 
single Polygon or Line objects, see the thread at:

https://stat.ethz.ch/pipermail/r-sig-geo/2009-February/005037.html

and my reply yesterday on stepping through Polygons and Polygon lists 
using lapply, and where to look for sample code. Basically:

all_the_coords <- lapply(slot(scot_BNG, "polygons"), function(Plns)
   lapply(slot(Plns, "Polygons"), function(Pln) slot(Pln, "coords")))

is a list of lists of coordinate matrices.

Hope this helps,

Roger

>
> Thanks!
>
> --j 
> Hengl, T. wrote:
>> 
>> If you work with large shapes/grids, try also SAGA:
>> 
>> > rsaga.get.usage(lib="grid_gridding", 3)
>> SAGA CMD 2.0.3
>> library path:   C:/Progra~1/saga_vc/modules
>> library name:   grid_gridding
>> module name :   Shapes to Grid
>> ...
>> 
>> But before SAGA, you need to reproject the polygons first (if necessary).
>> 
>> see also some examples from:
>> http://spatial-analyst.net/wiki/index.php?title=Export_maps_to_GE#Polygon_maps
>> 
>> HTH
>> 
>> Tom Hengl
>> 
>> -----Original Message-----
>> From: r-sig-geo-bounces at stat.math.ethz.ch on behalf of Paul Hiemstra
>> Sent: Sat 2/14/2009 10:07 PM
>> To: Jonathan Greenberg
>> Cc: r-sig-geo at stat.math.ethz.ch
>> Subject: Re: [R-sig-Geo] Vector to raster?
>> 
>> Jonathan Greenberg schreef:
>> > How do I take a polygon in some OGR supported vector layer (say, a
>> > shapefile), and rasterize it given a pixel size and projection/datum?
>> >
>> > --j
>> >
>> Take a look at the spsample() function.
>> 
>> Paul
>> 
>> --
>> Drs. Paul Hiemstra
>> Department of Physical Geography
>> Faculty of Geosciences
>> University of Utrecht
>> Heidelberglaan 2
>> P.O. Box 80.115
>> 3508 TC Utrecht
>> Phone:     +31302535773
>> Fax:    +31302531145
>> http://intamap.geo.uu.nl/~paul <http://intamap.geo.uu.nl/%7Epaul>
>> 
>> _______________________________________________
>> R-sig-Geo mailing list
>> R-sig-Geo at stat.math.ethz.ch
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>> 
>
>

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