[R-sig-Geo] Using the Z value for 3-D polygons?
Torleif Markussen Lunde
torleif.lunde at cih.uib.no
Thu Mar 5 16:51:20 CET 2009
Ok. Again I will not answer the question exactly the way you would like to
have things done. Anyway this is an approach to use lattice.
1. Use spsample to create a grid
2. Make a data.frame and give the corresponding data values from the shp-file
to the coordinates in the data frame
3. plot with wireframe from lattice.
4. An alternative approach
require(maptools)
require(lattice)
nc <- readShapePoly(system.file("shapes/sids.shp", package="maptools")[1],
proj4string=CRS("+proj=longlat +datum=NAD27"))
#Make a grid. Details controlled by cellsize
grd <- spsample(nc, cellsize = c(.05,.05),
type = "regular",
offset = c(0.5, 0.5))
#Get polygon ID and give "value" to z data
z <- overlay(as(grd, "SpatialPoints"), nc)
grd.co <- data.frame(slot(grd, 'coords'))
grd.df <- data.frame(cbind(grd.co, z))
grd.df$z <- sapply(grd.df$z,
function(i) slot(nc, 'data')$NWBIR79[i])
#Make the plot. Notice that wireframe is not drawn at NA values
with(grd.df, wireframe(z ~ x1 * x2,
scales = list(arrows = FALSE),
shade =TRUE, pretty = TRUE,
xlab="", ylab="", zlab=""))
#Trying to correct by creating matrix of z data
grd.sp <- SpatialPixelsDataFrame(cbind(grd.df$x1, grd.df$x2),
data = data.frame(grd.df$z))
xyz <- as.image.SpatialGridDataFrame(grd.sp["grd.df.z"])
xyz$z[is.na(xyz$z)] <- 0
with(xyz, wireframe(z, row.values=x, col.values=y, shade = TRUE,
light.source = c(90,190,90)),
scales = list(arrows = FALSE))
At least this is one approach. Since I am not too familiar with spplot, I
apologize if there should be an option to plot the polygons directly with z
in the third dimension. From what Edzer writes I would guess it is not
straight forward. In that case I would be interested to know.
The rendering is still not too pretty.
Best wishes
Torleif
On Thursday 05 March 2009 02:48:53 am Jim Burke wrote:
> Hi, I still need some help. Scant progress.
>
> Resources: I have Deepayan's "Lattice" book and Roger,
> Edzer, and Virgillio's "Applied Spatial Data Analysis"
> if you want to refer me to anything there.
>
> A. I use spplot with spplot.polygon parameters I receive the
> following error.
> . The argument seems there and returns "iso".
> . Am I missing something? My call is in section I below.
> "Error in append(list(formula, data = as(sdf, "data.frame"),
> aspect = aspect,: argument is missing, with no default"
>
> B. By the way, mapasp doc could be cleaner. In the line below...
> 1) no opening brace to match "180)" closing brace
> 2) "middle of the map" the correct word is probably "median" not "mean"
> "let s = 1/cos(My * pi)/180) with My the y coordinate of the middle
> of the map (the mean of ylim)"
>
> C. I would be happy to email my spatial polygon data frame if you want
> to experiment.
>
> Below is my non working spplot call. Idea is to have a 3d version
> that uses the z value.
>
> I. Errors with aspect function argument
> library(sp)
> library(maptools)
> library(rgdal)
> library(PBSmapping)
> tx3_sp <- readShapePoly("tx3_sp.shp", IDvar="PCT",
> proj4string=CRS("+proj=aea +ellps=GRS80 +datum=WGS84"))
>
> spplot (
> tx3_sp, # obj
> c("D_CNT"), # zcol = names(obj). D_CNT has a max of
> 1693 , # other arguments passed to levelplot (grids,
> # polygons) or xyplot (points)
> , # names, to use in panel, if different
> from zcol
> # names
> scales = list(draw = FALSE), # scales argument to be passed to
> Lattice
> # plots;
> xlab = NULL,
> ylab = NULL,
> aspect = mapasp(obj,xlim,ylim),
> # aspect ratio for spatial axes;
> defaults to "iso"
> # (one unit on the x-axis equals one
> unit on the
> # y-axis) but may be set to more
> suitable values
> # if the data are e.g. if coordinates are
> #latitude/longitude
> panel = panel.polygonsplot,
> # depending on the class of obj,
> panel.polygonsplot
> # (for polygons or lines),
> panel.gridplot (grids) or
> # panel.pointsplot (points) is used;
> for further
> # control custom panel functions
> # can be supplied that call one of
> these panel functions,
> # but do read how the argument
> sp.layout may help
> sp.layout = NULL, # NULL or list; see notes below
> , # formula, optional; may be useful to
> plot a transformed
> # value. Defaults to z~x+y for single
> and z~x+y|name
> # for multiple attributes; use e.g.
> exp(x)~x+y|name to
> # plot the exponent of the
> z-variable
> xlim = bbox(tx3_sp)[1, ],# numeric; x-axis limits
> ylim = bbox(tx3_sp)[2, ] # numeric; y-axis limits
> )
>
> II. > summary(tx3_sp)
> Object of class SpatialPolygonsDataFrame
> Coordinates:
> min max
> r1 2488376 2542115
> r2 7003736 7046215
> Is projected: TRUE
> proj4string :
> [+proj=aea +ellps=GRS80 +datum=WGS84 +lat_1=29.5 +lat_2=45.5
> +towgs84=0,0,0]
> Data attributes: (I did not list these)
>
> III. mapasp(tx3_sp,bbox(tx3_sp)[1, ],bbox(tx3_sp)[2, ])
> [1] "iso"
>
>
> Thanks,
> Jim Burke
>
> Jim Burke wrote:
> > Hi Everyone,
> >
> > QUESTION: How can I plot a spatial polygon (SpatialPolygonsDataFrame)
> > as 3-d? Perhaps using "spplot.polygons" with its z value?
> >
> > GOAL: To have a series of polygons look like approximately like the
> > last example in the PDF below.
> > http://ocw.mit.edu/NR/rdonlyres/Urban-Studies-and-Planning/11-521Spatial-
> >Database-Management-and-Advanced-Geographic-Information-SystemsSpring2003/
> >ACA80C9F-4089-403E-9669-50763381D08A/0/lect13c.pdf
> >
> >
> > EXAMPLE SUGGESTION: Use "spplot.pologons" with the venerable
> > SourceForge NC fig21.R example and using its color coding values as
> > the polygon $z value? Assuming lat and long are in the sp.
> > http://r-spatial.sourceforge.net/gallery/#fig21.R
> >
> > Then if its not as straightforward as above, what would be the steps
> > to convert a spatial polygon to points or something to graph that way?
> >
> > Thanks,
> > Jim Burke
> >
> > _______________________________________________
> > R-sig-Geo mailing list
> > R-sig-Geo at stat.math.ethz.ch
> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
More information about the R-sig-Geo
mailing list