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-boun...@stat.math.ethz.ch on behalf of Paul Hiemstra
Sent: Sat 2/14/2009 10:07 PM
To: Jonathan Greenberg
Cc: r-sig-geo@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@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: jgrn...@hotmail.com, Gchat: jgrn307
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
