[R-sig-Geo] Vector to raster?

[Jonathan Greenberg]

Roger:

    Thanks, I'll try to incorporate some of your suggestions and post 
the code + some sample imagery next time I ask a question.  I'm trying 
to build up some R-based code that replicates some of what my "starspan" 
program does (starspan.casil.ucdavis.edu) -- essentially taking vector 
data, querying out the raster data at those locations (point, poly or 
line), and fusing the vector DB with the extracted raster data.  We'll 
probably eventually work on the C version of the code again, and write 
an R wrapper, but for now I wanted an R-only script.

    The larger idea is to be able to generate training/test data for use 
in remote sensing classification and analysis.  We're focusing on using 
rgdal for a course at UC Davis on open source geospatial analysis.  I 
just got your book, so I'll be reading it this weekend!

--j

Roger Bivand wrote:
> 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
>>>
>>
>>
>

-- 

Jonathan A. Greenberg, PhD
Postdoctoral Scholar
Center for Spatial Technologies and Remote Sensing (CSTARS)
University of California, Davis
One Shields Avenue
The Barn, Room 250N
Davis, CA 95616
Cell: 415-794-5043
AIM: jgrn307, MSN: jgrn307 at hotmail.com, Gchat: jgrn307