[R-sig-Geo] library(shapefiles): write PolyLineZ or PolygonZ shapes?

[Albrecht Gebhardt]

hi,

I'm trying to write some 3D polygoins as ESRI shapefile.

While write.shp() has some code for type==13 or 15 (which is PolyLineZ
and PolygonZ), there is no such part in convert.to.shapefile().

example:

first generate some X,Y,Z data and triangulate it (only X and Y):

library(tripack)
xyz<-as.data.frame(cbind(runif(5),runif(5),rnorm(5)))
colnames(xyz)<-c("x","y","z")
t<-tri.mesh(xyz$x,xyz$y)
tr<-triangles(t)
# generate one polyline per triangle:
Id <- NULL; X <- NULL; Y <- NULL; Z <- NULL
for(i in 1:dim(tr)[1]) {
  Id <- c(Id,rep(i,3))
  X <- c(X, xyz$x[tr[i,1]], xyz$x[tr[i,2]], xyz$x[tr[i,3]])
  Y <- c(Y, xyz$y[tr[i,1]], xyz$y[tr[i,2]], xyz$y[tr[i,3]])
  Z <- c(Z, xyz$z[tr[i,1]], xyz$z[tr[i,2]], xyz$z[tr[i,3]])
}
dd <- data.frame(Id=Id,X=X,Y=Y, Z=Z)

dd contains now something like (three triangles):

  Id          X         Y          Z
1  1 0.47578460 0.7062264  1.7720352
2  1 0.07356501 0.6989611 -0.3850338
3  1 0.24476600 0.2560522  0.5450287
4  2 0.47578460 0.7062264  1.7720352
5  2 0.24476600 0.2560522  0.5450287
6  2 0.47912981 0.2355883 -0.2675531
7  3 0.47578460 0.7062264  1.7720352
8  3 0.47912981 0.2355883 -0.2675531
9  3 0.78672806 0.3244467  0.8620736

If I now try (in analogy to the write.shapefile PolyLine example)

library(shapefiles)
ddTable <- data.frame(Id=1:dim(tr)[1], Name=as.character(1:dim(tr)[1]))
ddShapefile <- convert.to.shapefile(dd, ddTable, "Id", 13)
write.shapefile(ddShapefile, "~/test/tri", arcgis=TRUE)

I get an error:

Fehler in "dimnames<-.data.frame"(`*tmp*`, value = list(c("1", "2", "3",  :
    ungültige 'dimnames' für data frame gegeben   # invalid dimnames

Now I tried to modify convert.to.shapefile() (see attachement) so that
it somehow handles the type==13/15 case (but it contains some blind
guesses regarding the content.sizes, I should read
http://www.esri.com/library/whitepapers/pdfs/shapefile.pdf more
precisely) and then at least it runs:

ddShapefile <- my.convert.to.shapefile(dd, ddTable, "Id", 13)
write.shapefile(ddShapefile, "~/test/tri", arcgis=TRUE)

I can even re-read this file with

f <- read.shapefile("~/test/tri")

and "f" looks identical (but I didnt compare it bit by bit) to
"ddShapefile"

But GRASS doesn't like it, it only reads the first triangle:

GRASS 6.0.1 (import_xy):~ > v.in.ogr dsn=~/test/tri.shp
output=tri_import min_area=0.0001 snap=-1 -z
Projection of input dataset and current location appear to match.
Proceeding with import...
Layer: tri
WARNING: 2 features without geometry.
-----------------------------------------------------
Building topology ...
1 primitives registered
Building areas:  100%
0 areas built
0 isles built
Attaching islands:
Attaching centroids:  100%
Topology was built.
Number of nodes     :   2
Number of primitives:   1
Number of points    :   0
Number of lines     :   1
Number of boundaries:   0
Number of centroids :   0
Number of areas     :   0
Number of isles     :   0





So my question is:

Is it really necessary to extend functionality of convert.to.shapefile()
to be able to write 3D polygons/polylines or am I simply missing the
correct function?



Thanks

Albrecht
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: my.convert.to.shapefile.R
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20070710/02383d65/attachment.pl>