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

[Roger Bivand]

On Tue, 10 Jul 2007, Albrecht Gebhardt wrote:

> hi,
>
> I'm trying to write some 3D polygons as ESRI shapefile.

I think that this is quite hazardous. We do see 3D shapefiles in the wild, 
usually with Z=0 for all coordinates. I guess that something could be done 
in maptools, but not using sp classes, which are 2D for lines and 
polygons. I'm not sure if the shapefiles package author tracks this list 
regularly, so I've added him to this reply.

Roger

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

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