Nothing
R/shapefiles.R
In shapefiles: Read and Write ESRI Shapefiles
Defines functions
dp
change.id
convert.to.simple
convert.to.shapefile
scaleXY
add.xy
write.shapefile
read.shapefile
write.dbf
read.dbf
write.shx
calc.header
read.shx
write.shp
read.shp
Documented in
add.xy
calc.header
change.id
convert.to.shapefile
convert.to.simple
dp
read.dbf
read.shapefile
read.shp
read.shx
scaleXY
write.dbf
write.shapefile
write.shp
write.shx
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Shapefile Format - Read/Write shapefile format within R
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Created 2/13/03 Ben Stabler benjamin.stabler@odot.state.or.us
# Revised 6/11/03 Ben Stabler benjamin.stabler@odot.state.or.us
# Revised 7/7/03 Ben Stabler benjamin.stabler@odot.state.or.us
# Revised 7/23/03 Ben Stabler benjamin.stabler@odot.state.or.us
# Revised 8/10/05 Ben Stabler, benstabler@yahoo.com
# Reading of shape type 13 and 15 from Don MacQueen, macq@llnl.gov
# Conversion of simple polygons to shapefile format
# from Manuel Chirouze, Manuel.Chirouze@benfieldgroup.com
# Small bug fix to write.dbf to deal with 1 column white space correctly
# Revised 8/15/05 Ben Stabler, to use write.dbf and read.dbf from foreign library
# Added dp function as well
# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU General Public License
# as published by the Free Software Foundation; either version 2
# of the License, or (at your option) any later version.
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# DESCRIPTION
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#source("shapefiles.R")
#myshape <- read.shapefile("myshapewithnoextension")
#myshape$shp - the shp file
#myshape$shx - the index file
#myshape$dbf - the dbf file
#Each of the three lists of the myshape contain a header and some data. So there is...
#myshape$shp$header (the header info) and myshape$shp$shp (the geographic data)
#myshape$shx$header and myshape$shx$index (the index data)
#myshape$dbf$header and myshape$dbf$dbf (the dbf stored as a R data.frame)
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#Read in a point, line or polygon SHP file
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
read.shp <- function(shp.name) {
infile<-file(shp.name,"rb")
#Header
file.code <- readBin(infile,integer(),1, endian="big")
unused <- readBin(infile,integer(),5, endian="big")
file.length <- readBin(infile,integer(),1, endian="big")
file.version <- readBin(infile,integer(),1, endian="little")
shape.type <- readBin(infile,integer(),1, endian="little")
xmin <- readBin(infile,double(),1, endian="little")
ymin <- readBin(infile,double(),1, endian="little")
xmax <- readBin(infile,double(),1, endian="little")
ymax <- readBin(infile,double(),1, endian="little")
zmin <- readBin(infile,double(),1, endian="little")
zmax <- readBin(infile,double(),1, endian="little")
mmin <- readBin(infile,double(),1, endian="little")
mmax <- readBin(infile,double(),1, endian="little")
header <- list(file.code,file.length,file.version,shape.type,xmin,ymin,xmax,ymax,zmin,zmax,mmin,mmax)
names(header) <- c("file.code","file.length","file.version","shape.type","xmin","ymin","xmax","ymax","zmin","zmax","mmin","mmax")
rm(file.code,file.length,file.version,shape.type,xmin,ymin,xmax,ymax,zmin,zmax,mmin,mmax,unused)
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#Point Shapefile
if (header$shape.type==1) {
record <- integer()
shape.type <- integer()
x <- double()
y <- double()
record.num <- 1
while (length(record.num)!=0) {
#Record Number
record.num <- readBin(infile,integer(),1, 4, endian="big")
if (length(record.num)==0) break
record <- c(record,record.num)
#Content Length
content.length <- readBin(infile,integer(),1, 4, endian="big")
#Shape Type (should be 1)
temp.type <- readBin(infile,integer(),1, 4, endian="little")
shape.type <- c(shape.type, temp.type)
if (temp.type==1) {
#X and Y values
x <- c(x,readBin(infile,double(),1, 8, endian="little"))
y <- c(y,readBin(infile,double(),1, 8, endian="little"))
#To read in extra data created by GeoMedia shapefile export
if (content.length>10) {
temp <- readBin(infile,integer(), ((content.length-10)/2), 4, endian="big")
}
}
if (temp.type==0) {
#Null shape type has no data
x <- c(x,NA)
y <- c(y,NA)
}
}
point.data <- cbind(record=record, x=x, y=y, shape.type=shape.type)
close(infile)
}
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#PolyLine, Polygon, PolyLineZ, and PolygonZ Shapefile
if (header$shape.type==3 || header$shape.type==5 || header$shape.type == 13 || header$shape.type == 15) {
shape <- list()
record.num <- 1
while (length(record.num)!=0) {
#Record Number
record.num <- readBin(infile,integer(),1, endian="big")
if (length(record.num)==0) break
#Content Length
content.length <- readBin(infile,integer(),1, endian="big")
#Shape Type
shape.type <- readBin(infile, integer(), 1, endian="little")
if (shape.type==3 || shape.type==5) {
box <- readBin(infile, double(), 4, endian="little")
names(box) <- c("xmin","ymin","xmax","ymax")
num.parts <- readBin(infile, integer(), 1, endian="little")
num.points <- readBin(infile, integer(), 1, endian="little")
parts <- readBin(infile, integer(), num.parts, endian="little")
points <- readBin(infile, double(), num.points*2, endian="little")
X <- points[seq(1,(num.points*2),by=2)]
Y <- points[seq(2,(num.points*2),by=2)]
points <- data.frame(X,Y)
}
#Type 13 polylineZ (and 15 polygonZ since it is the same as 13)
if (shape.type==13 || shape.type==15) {
box <- readBin(infile, double(), 4, endian = "little")
names(box) <- c("xmin", "ymin", "xmax", "ymax")
num.parts <- readBin(infile, integer(), 1, endian = "little")
num.points <- readBin(infile, integer(), 1, endian = "little")
parts <- readBin(infile, integer(), num.parts, endian = "little")
points <- readBin(infile, double(), num.points * 2, endian = "little")
X <- points[seq(1, (num.points * 2), by = 2)]
Y <- points[seq(2, (num.points * 2), by = 2)]
#Additional information for type 13 and 15
zrange <- readBin(infile,double(),2,endian='little')
Z <- readBin(infile,double(),num.points,endian='little')
mrange <- readBin(infile,double(),2,endian='little')
M <- readBin(infile,double(),num.points,endian='little')
points <- data.frame(X, Y, Z, M)
}
if (shape.type==0) {
#Null shape type has no data
box <- rep(NA,4)
num.parts <- NA
num.points <- NA
parts <- NA
points <- data.frame(X=NA,Y=NA)
}
shape.info <- list(record=record.num, content.length=content.length, shape.type=shape.type, box=box, num.parts=num.parts, num.points=num.points, parts=parts, points=points)
shape <- c(shape,list(shape.info))
}
close(infile)
}
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if (header$shape.type==1) {
return.val <- point.data
}
if (header$shape.type==3 || header$shape.type==5 || header$shape.type == 13 || header$shape.type == 15) {
return.val <- shape
}
list(shp=return.val,header=header)
}
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#Write out a point, line or polygon SHP file
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
write.shp <- function(shp, out.name) {
outfile<-file(out.name,"wb")
#Header
writeBin(as.integer(9994), outfile, endian="big")
writeBin(as.integer(rep(0,5)), outfile, endian="big")
writeBin(as.integer(shp$header$file.length), outfile, endian="big")
writeBin(as.integer(1000), outfile, endian="little")
writeBin(as.integer(shp$header$shape.type), outfile, endian="little")
writeBin(as.double(shp$header$xmin), outfile, endian="little")
writeBin(as.double(shp$header$ymin), outfile, endian="little")
writeBin(as.double(shp$header$xmax), outfile, endian="little")
writeBin(as.double(shp$header$ymax), outfile, endian="little")
writeBin(as.double(shp$header$zmin), outfile, endian="little")
writeBin(as.double(shp$header$zmax), outfile, endian="little")
writeBin(as.double(shp$header$mmin), outfile, endian="little")
writeBin(as.double(shp$header$mmax), outfile, endian="little")
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if (shp$header$shape.type==1) {
for (record in 1:nrow(shp$shp)) {
if (shp$shp[record,4]==1) {
#Record Number
writeBin(as.integer(shp$shp[record,1]), outfile, endian="big")
#Content Length
writeBin(as.integer(10), outfile, endian="big")
#Shape Type
writeBin(as.integer(1), outfile, endian="little")
#X - For Points only
writeBin(as.double(shp$shp[record,2]), outfile, endian="little")
#Y - For Points only
writeBin(as.double(shp$shp[record,3]), outfile, endian="little")
}
if (shp$shp[record,4]==0) {
#Record Number
writeBin(as.integer(shp$shp[record,1]), outfile, endian="big")
#Content Length
writeBin(as.integer(4), outfile, endian="big")
#Shape Type
writeBin(as.integer(0), outfile, endian="little")
}
}
}
if (shp$header$shape.type==3 || shp$header$shape.type==5 || shp$header$shape.type == 13 || shp$header$shape.type == 15) {
for (record in 1:length(shp$shp)) {
if (shp$shp[[record]]$shape.type==3 || shp$shp[[record]]$shape.type==5) {
#Record Number
writeBin(as.integer(shp$shp[[record]]$record), outfile, endian="big")
#Content Length
writeBin(as.integer(shp$shp[[record]]$content.length), outfile, endian="big")
#Shape Type
writeBin(as.integer(shp$shp[[record]]$shape.type), outfile, endian="little")
#Bounding Box
writeBin(as.double(shp$shp[[record]]$box[1]), outfile, endian="little")
writeBin(as.double(shp$shp[[record]]$box[2]), outfile, endian="little")
writeBin(as.double(shp$shp[[record]]$box[3]), outfile, endian="little")
writeBin(as.double(shp$shp[[record]]$box[4]), outfile, endian="little")
#Number of parts (segments)
writeBin(as.integer(shp$shp[[record]]$num.parts), outfile, endian="little")
#Number of total points
writeBin(as.integer(shp$shp[[record]]$num.points), outfile, endian="little")
#Parts
writeBin(as.integer(shp$shp[[record]]$parts), outfile, endian="little")
#Points - merge the X and Y points into one vector
point.stream <- rep(0,length(shp$shp[[record]]$points$X)*2)
point.stream[seq(1,length(shp$shp[[record]]$points$X)*2,by=2)] <- shp$shp[[record]]$points$X
point.stream[seq(2,length(shp$shp[[record]]$points$Y)*2,by=2)] <- shp$shp[[record]]$points$Y
writeBin(as.double(point.stream), outfile, endian="little")
#Need to check to make sure first and last vertex are the same for Polygon
}
#Type 13 polylineZ (and 15 polygonZ since it is the same as 13)
if (shp$shp[[record]]$shape.type==13 || shp$shp[[record]]$shape.type==15) {
#Record Number
writeBin(as.integer(shp$shp[[record]]$record), outfile, endian="big")
#Content Length
writeBin(as.integer(shp$shp[[record]]$content.length), outfile, endian="big")
#Shape Type
writeBin(as.integer(shp$shp[[record]]$shape.type), outfile, endian="little")
#Bounding Box
writeBin(as.double(shp$shp[[record]]$box[1]), outfile, endian="little")
writeBin(as.double(shp$shp[[record]]$box[2]), outfile, endian="little")
writeBin(as.double(shp$shp[[record]]$box[3]), outfile, endian="little")
writeBin(as.double(shp$shp[[record]]$box[4]), outfile, endian="little")
#Number of parts (segments)
writeBin(as.integer(shp$shp[[record]]$num.parts), outfile, endian="little")
#Number of total points
writeBin(as.integer(shp$shp[[record]]$num.points), outfile, endian="little")
#Parts
writeBin(as.integer(shp$shp[[record]]$parts), outfile, endian="little")
#Points - merge the X and Y points into one vector
point.stream <- rep(0,length(shp$shp[[record]]$points$X)*2)
point.stream[seq(1,length(shp$shp[[record]]$points$X)*2,by=2)] <- shp$shp[[record]]$points$X
point.stream[seq(2,length(shp$shp[[record]]$points$Y)*2,by=2)] <- shp$shp[[record]]$points$Y
writeBin(as.double(point.stream), outfile, endian="little")
#Need to check to make sure first and last vertex are the same for Polygon
#Write out Z Range and Z values
minZ <- min(shp$shp[[record]]$points$Z)
maxZ <- max(shp$shp[[record]]$points$Z)
writeBin(as.double(minZ), outfile, endian='little')
writeBin(as.double(maxZ), outfile, endian='little')
writeBin(as.double(shp$shp[[record]]$points$Z), outfile, endian='little')
#Write out M Range and M values
minM <- min(shp$shp[[record]]$points$M)
maxM <- max(shp$shp[[record]]$points$M)
writeBin(as.double(minM), outfile, endian='little')
writeBin(as.double(maxM), outfile, endian='little')
writeBin(as.double(shp$shp[[record]]$points$M), outfile, endian='little')
}
if (shp$shp[[record]]$shape.type==0) {
#Record Number
writeBin(as.integer(shp$shp[[record]]$record), outfile, endian="big")
#Content Length
writeBin(as.integer(shp$shp[[record]]$content.length), outfile, endian="big")
#Shape Type
writeBin(as.integer(shp$shp[[record]]$shape.type), outfile, endian="little")
}
}
}
close(outfile)
}
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#Read in the SHX file (index)
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
read.shx <- function(shx.name) {
infile<-file(shx.name,"rb")
#Header
file.code <- readBin(infile,integer(),1, endian="big")
unused <- readBin(infile,integer(),5, endian="big")
file.length <- readBin(infile,integer(),1, endian="big")
file.version <- readBin(infile,integer(),1, endian="little")
shape.type <- readBin(infile,integer(),1, endian="little")
xmin <- readBin(infile,double(),1, endian="little")
ymin <- readBin(infile,double(),1, endian="little")
xmax <- readBin(infile,double(),1, endian="little")
ymax <- readBin(infile,double(),1, endian="little")
zmin <- readBin(infile,double(),1, endian="little")
zmax <- readBin(infile,double(),1, endian="little")
mmin <- readBin(infile,double(),1, endian="little")
mmax <- readBin(infile,double(),1, endian="little")
header <- list(file.code,file.length,file.version,shape.type,xmin,ymin,xmax,ymax,zmin,zmax,mmin,mmax)
names(header) <- c("file.code","file.length","file.version","shape.type","xmin","ymin","xmax","ymax","zmin","zmax","mmin","mmax")
records <- (file.length-50)/4
#Record Data
index.data <- matrix(0,records,2)
colnames(index.data) <- c("Offset","Length")
for (record in 1:records) {
#First record offset is 50
index.data[record,1] <- readBin(infile,integer(),1, endian="big")
#Content length
index.data[record,2] <- readBin(infile,integer(),1, endian="big")
}
close(infile)
list(index=index.data,header=header)
}
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Function to calculate header information for a point file
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Used to clean up GeoMedia shapefile export files
calc.header <- function(shapefile) {
#1) Calculate shp header info
#Calculate number of bytes for null record headers
num.null.record.bytes <- length(shapefile$shp$shp[,2][is.na(shapefile$shp$shp[,2])]) * 4
#Calculate number of bytes for complete record record headers
num.record.header.bytes <- length(shapefile$shp$shp[,2]) * 8
#Calculate number of bytes for records (both geoMedia and generic shapefile)
#Should be value 20 (56 for GeoMedia files)
num.record.bytes <- (length(shapefile$shp$shp[,2]) - length(shapefile$shp$shp[,2][is.na(shapefile$shp$shp[,2])])) * 20
#Sum the bytes to get total bytes (100 = header bytes)
total.bytes <- 100 + sum(num.null.record.bytes, num.record.header.bytes, num.record.bytes)
#Divide by 2 to get number of 16-bit words
file.length <- total.bytes / 2
#Replace file.length in shapefile header with new file.length calculation
shapefile$shp$header$file.length <- file.length
#Replace lengths of 28 with 10 for shx from GeoMedia
#2) Calculate shx info
if (shapefile$shx$index[1,2]==28) {
shapefile$shx$index[,2][shapefile$shx$index[,2]==28] <- 10
}
#Set the content.length to 14 in order to calculate offset
shapefile$shx$index[,2][shapefile$shx$index[,2]==10]<-14
#Set null value content.lengths from 2 to 6 to calculate offsets
shapefile$shx$index[,2][shapefile$shx$index[,2]==2]<-6
#Calculate new offset values for each record from beginning of shp file
new.offset <- cumsum(c(50,shapefile$shx$index[,2]))
#Replace the existing offsets with the new ones
shapefile$shx$index[,1] <- new.offset[1:(length(new.offset)-1)]
#Set the content.lenghts back to their respective values
shapefile$shx$index[,2][shapefile$shx$index[,2]==14]<-10
shapefile$shx$index[,2][shapefile$shx$index[,2]==6]<-2
shapefile
}
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#Write out SHX file (index)
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
write.shx <- function(shx, out.name) {
outfile<-file(out.name,"wb")
#Header
writeBin(as.integer(9994), outfile, endian="big")
writeBin(as.integer(rep(0,5)), outfile, endian="big")
writeBin(as.integer(shx$header$file.length), outfile, endian="big")
writeBin(as.integer(1000), outfile, endian="little")
writeBin(as.integer(shx$header$shape.type), outfile, endian="little")
writeBin(as.double(shx$header$xmin), outfile, endian="little")
writeBin(as.double(shx$header$ymin), outfile, endian="little")
writeBin(as.double(shx$header$xmax), outfile, endian="little")
writeBin(as.double(shx$header$ymax), outfile, endian="little")
writeBin(as.double(shx$header$zmin), outfile, endian="little")
writeBin(as.double(shx$header$zmax), outfile, endian="little")
writeBin(as.double(shx$header$mmin), outfile, endian="little")
writeBin(as.double(shx$header$mmax), outfile, endian="little")
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
for (record in 1:nrow(shx$index)) {
#Record Offset
writeBin(as.integer(shx$index[record,1]), outfile, endian="big")
#Content Length
writeBin(as.integer(shx$index[record,2]), outfile, endian="big")
}
close(outfile)
}
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#Read DBF format
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
read.dbf <- function(dbf.name, header=FALSE) {
#Load foreign package for read.dbf
## library(foreign)
#Read header if needed
if(header) {
infile<-file(dbf.name,"rb")
#Header
file.version <- readBin(infile,integer(), 1, size=1, endian="little")
file.year <- readBin(infile,integer(), 1, size=1, endian="little")
file.month <- readBin(infile,integer(), 1, size=1, endian="little")
file.day <- readBin(infile,integer(), 1, size=1, endian="little")
num.records <- readBin(infile,integer(), 1, size=4, endian="little")
header.length <- readBin(infile,integer(), 1, size=2, endian="little")
record.length <- readBin(infile,integer(), 1, size=2, endian="little")
file.temp <- readBin(infile,integer(), 20, size=1, endian="little")
header <- list(file.version,file.year, file.month, file.day, num.records, header.length, record.length)
names(header) <- c("file.version","file.year","file.month","file.day","num.records","header.length","record.length")
rm(file.version,file.year, file.month, file.day, num.records, header.length, record.length)
#Calculate the number of fields
num.fields <- (header$header.length-32-1)/32
field.name <- NULL
field.type <- NULL
field.length <- NULL
field.decimal <- NULL
#Field Descriptions (32 bytes each)
for (i in 1:num.fields) {
field.name.test <- readBin(infile,character(), 1, size=10, endian="little")
field.name <- c(field.name,field.name.test)
if (nchar(field.name.test)!=10) {
file.temp <- readBin(infile,integer(), 10-(nchar(field.name.test)), 1, endian="little")
}
field.type <- c(field.type,readChar(infile, 1))
file.temp <- readBin(infile,integer(), 4, 1, endian="little")
field.length <- c(field.length,readBin(infile,integer(), 1, 1, endian="little"))
field.decimal <- c(field.decimal, readBin(infile,integer(), 1, 1, endian="little"))
file.temp <- readBin(infile,integer(), 14, 1, endian="little")
}
#Create a table of the field info
fields <- data.frame(NAME=field.name,TYPE=field.type,LENGTH=field.length,DECIMAL=field.decimal)
#Set all fields with length<0 equal to correct number of characters
fields$LENGTH[fields$LENGTH<0]<-(256+fields$LENGTH[fields$LENGTH<0])
#Read in end of attribute descriptions terminator - should be integer value 13
file.temp <- readBin(infile,integer(), 1, 1, endian="little")
#Increase the length of field 1 by one to account for the space at the beginning of each record
fields$LENGTH[1]<-fields$LENGTH[1]+1
#Add fields to the header list
header <- c(header,fields=NULL)
header$fields <- fields
#Close file connection
close(infile)
}
#Read in dbf
dbf <- foreign::read.dbf(dbf.name)
#Return the dbf as a list with a data.frame and a header list
list(dbf=dbf, header=header)
}
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#Write out DBF format
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
write.dbf <- function(dbf, out.name, arcgis=FALSE) {
#Load foreign package for write.dbf
## library(foreign)
#If output intended for ArcGIS, replace "." with "_" for column names in header
if (arcgis==T) {
colnames(dbf$dbf) <- gsub("\\.","_",colnames(dbf$dbf))
}
#Write out dbf
foreign::write.dbf(dbf$dbf, out.name)
}
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#Function to Read in a Shapefile
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
read.shapefile <- function(shape.name) {
shp.data <- read.shp(paste(shape.name, ".shp", sep=""))
shx.data <- read.shx(paste(shape.name, ".shx", sep=""))
dbf.data <- get("read.dbf","package:shapefiles")(paste(shape.name, ".dbf", sep=""))
shapefile <- list(shp=shp.data,shx=shx.data,dbf=dbf.data)
shapefile
}
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#Function to Write out a Shapefile
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
write.shapefile <- function(shapefile, out.name, arcgis=FALSE) {
write.shp(shapefile$shp, paste(out.name, ".shp", sep=""))
write.shx(shapefile$shx, paste(out.name, ".shx", sep=""))
get("write.dbf","package:shapefiles")(shapefile$dbf, paste(out.name, ".dbf", sep=""), arcgis)
}
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#Function to Add the X and Y Coordinates to the DBF of a shapefile
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
add.xy <- function(shapefile) {
if (shapefile$shp$header$shape.type !=1) {
stop ("Must be point shapefile")
}
shapefile$dbf[[1]]$XCOORD <- shapefile$shp[[1]][,2]
shapefile$dbf[[1]]$YCOORD <- shapefile$shp[[1]][,3]
shapefile
}
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#Function to scale the X and Y Coordinates of the shapefile
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
scaleXY <- function(shapefile, scale.factor) {
#Scale the bounding box of the shapefile by the scale factor
shapefile$shp$header$xmin <- shapefile$shp$header$xmin/scale.factor
shapefile$shp$header$ymin <- shapefile$shp$header$ymin/scale.factor
shapefile$shp$header$xmax <- shapefile$shp$header$xmax/scale.factor
shapefile$shp$header$ymax <- shapefile$shp$header$ymax/scale.factor
#Scale the bounding box of the shx file by the scale factor
shapefile$shx$header$xmin <- shapefile$shx$header$xmin/scale.factor
shapefile$shx$header$ymin <- shapefile$shx$header$ymin/scale.factor
shapefile$shx$header$xmax <- shapefile$shx$header$xmax/scale.factor
shapefile$shx$header$ymax <- shapefile$shx$header$ymax/scale.factor
#Scale the X and Y if a point shape
if (shapefile$shp$header$shape.type == 1) {
shapefile$shp[[1]][,2] <- shapefile$shp[[1]][,2]/scale.factor
shapefile$shp[[1]][,3] <- shapefile$shp[[1]][,3]/scale.factor
}
#Scale the X and Y point values and bounding box if a line or polyon shape
if (shapefile$shp$header$shape.type==3 || shapefile$shp$header$shape.type==5) {
for (shape in 1:length(shapefile$shp$shp)) {
shapefile$shp[[1]][[shape]]$points[,1] <- shapefile$shp[[1]][[shape]]$points[,1]/scale.factor
shapefile$shp[[1]][[shape]]$points[,2] <- shapefile$shp[[1]][[shape]]$points[,2]/scale.factor
shapefile$shp[[1]][[shape]]$box <- shapefile$shp[[1]][[shape]]$box/scale.factor
}
}
shapefile
}
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#Converts a simple point, polyLine, or polygon data frame into a shapefile.
#The shpTable data frame must have three columns in this order: Id, X, and Y.
#The attTable data frame must have a field to join on - identified by field
#The type argument is either 1 (point), 3 (polyline) or 5 (polygon).
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
convert.to.shapefile <- function(shpTable, attTable, field, type) {
#Check to ensure field in attTable
if(!(field %in% colnames(attTable))) {
stop("Field not in attTable column names")
}
#Check to ensure the same Id set in shpTable and attTable
if(length(unique(shpTable[,1])) != length(attTable[,field])) {
stop("Different number of unique Ids in shpTable versus attTable")
}
if(any(sort(unique(shpTable[,1])) != sort(attTable[,field]))) {
stop("Id set in shpTable not the same as Id set in attTable")
}
#Build dbf - no longer uses header since "foreign.write.dbf" does not
dbf <- list(dbf=attTable,header=NULL)
#Build shp file
#shp list
shpTable[,c(2,3)] <- round(shpTable[,c(2,3)], 5)
colnames(shpTable) <- c("Id","X","Y")
v.content.length <- vector()
#point
if(type == 1) {
#loop for each entry in the dbf and build a shape
for (i in 1:nrow(attTable)) {
#Get shape entries for attTable Id i
shpTable.i <- shpTable[which(shpTable[,1] == attTable[i,field]), -1]
p.content.length <- 10 #shape.type and 2 doubles
v.content.length[i] <- p.content.length
if(i == 1) {
shp <- data.frame(
record = i,
x = shpTable.i[,1],
y = shpTable.i[,2],
shape.type = type)
} else {
shp2 <- data.frame(
record = i,
x = shpTable.i[,1],
y = shpTable.i[,2],
shape.type = type)
shp <- rbind(shp, shp2)
rm(shp2)
}
}
}
#polyLine and polygon shp
if(type == 3 || type == 5) {
#shp file
shp<-list()
#loop for each entry in the dbf and build a shape
for (i in 1:nrow(attTable)) {
#Get shape entries for attTable Id i
shpTable.i <- shpTable[which(shpTable[,1] == attTable[i,field]), -1]
num.points <- nrow(shpTable.i)
v.box <- c(min(shpTable.i[,1]),min(shpTable.i[,2]),max(shpTable.i[,1]),max(shpTable.i[,2]))
names(v.box) <- c("xmin","ymin","xmax","ymax")
p.content.length <- 24 + (num.points * 8)
v.content.length[i] <- p.content.length
shp[[i]]<-list(record = i,
content.length = p.content.length,
shape.type = type,
box = v.box,
num.parts = 1, #Only 1 part
num.points = nrow(shpTable.i),
parts = 0,
points = shpTable.i)
}
}
#Collect up all shp info
header<-list(file.code = 9994,
file.length = 50 + sum(v.content.length + 4),
file.version = 1000,
shape.type = type,
xmin = min(shpTable[,2]),
ymin = min(shpTable[,3]),
xmax = max(shpTable[,2]),
ymax = max(shpTable[,3]),
zmin = 0,
zmax = 0,
mmin = 0,
mmax = 0)
shp <- list(shp=shp, header=header)
#Build shx file
#shx list
shx<-list(index = cbind(Offset=c(50 ,(cumsum(v.content.length + 4) + 50)[-length(v.content.length)]),
Length = v.content.length),
header = shp$header)
shx$header$file.length <- 50 + 4 * nrow(shx$index)
#Return shapefile
shp.file <- list(shp=shp, shx=shx, dbf=dbf)
shp.file
}
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#Convert read.shp result to simple data.frame format (Id, X, Y)
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
convert.to.simple <- function(shp) {
#Determine type of shp
if(shp$header$shape.type == 1) {
#Get points
allPoints <- shp$shp[,c("record","x","y")]
colnames(allPoints) <- c("Id", "X", "Y")
} else {
#Get points
allPoints <- do.call(rbind, lapply(shp$shp, function(x) x$points[,c("X","Y")]))
#Calculate index and add to points
indexes <- 1:length(shp$shp)
repNum <- lapply(shp$shp, function(x) nrow(x$points))
allPoints <- cbind(Id=rep(indexes, repNum), allPoints)
}
#Return allPoints
allPoints
}
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#Change Index to a Field from the DBF
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#Changes the Index field of the shpTable to the new vector of fields
#Will replace the records in order, so Index 1 = element 1 in the new vector
change.id <- function(shpTable, newFieldAsVector) {
#Check to ensure the same Id set in shpTable and attTable
if(length(unique(shpTable[,1])) != length(newFieldAsVector)) {
stop("Different number of unique Ids in shpTable versus new field")
}
#Create new vector
count <- table(shpTable[,1])
shpTable[,1] <- rep(newFieldAsVector, count)
#Return result
shpTable
}
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#Douglas-Peucker Polyline Simplification
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#Ben Stabler, bstabler@ptvamerica.com, March 2005
#Douglas, D. and Peucker, T. (1973). "Algorithms for
#the reduction of the number of points required to
#represent a digitized line or its caricature."
#The Canadian Cartographer 10(2). 112-122.
#Currently uses the line, not the line segment to
#determine the distance of the points from the line.
#See http://www.lgc.com/resources/Doug_Peucker.pdf
#for more information. This can result in the
#omission of extreme "outlier-like" points. Try:
#points<-list(x=c(200,0,400),y=c(0,20,0))
#plot(points, type="l")
#lines(dp(points, 50), type="l", col="blue")
#Simple Example
#x <- c(5,3,4,1,8,9,10,11)
#y <- c(6,4,2,1,1,5,2,3)
#points <- list(x=x,y=y)
#plot(points, type="l")
#lines(dp(points, 2), type="l", col="blue")
####################################################
dp <- function(points, tolerance) {
#Convert to lowercase
names(points) <- tolower(names(points))
#Calculate distance between two points
distance <- function(x1,x2,y1,y2) {
sqrt((x2-x1)^2 + (y2-y1)^2)
}
#Calculate equation of a line from two points
equationOfLine <- function(x1,x2,y1,y2) {
slope <- (y2-y1)/(x2-x1)
b <- y1 - slope * x1
c(slope,b)
}
#Calculate y-intercept from a point and a slope
calcB <- function(slope, px, py) {
py - (slope * px)
}
#Calculate intercept of two lines from their slope and y-intercept
intercept <- function(s1,b1,s2,b2) {
x1 <- (b2-b1)/(s1-s2)
y1 <- s1 * x1 + b1
c(x1,y1)
}
#Setup vector to mark points to keep
keep <- rep(F, length(points$x))
keep[1] <- T
keep[length(keep)] <- T
#Function definition to simplify points
simplify <- function(start, end, tol=tolerance) {
#Calculate intermediate point distances
if (length(points$x[start:end]) > 2) {
#Avoid Inf slope
if( points$x[start] == points$x[end] ) { points$x[start] <- points$x[start] - 0.0000001 }
if( points$y[start] == points$y[end]) { points$y[start] <- points$y[start] - 0.0000001 }
#Calculate equation of line of middle points
line <- equationOfLine( points$x[start], points$x[end], points$y[start], points$y[end])
#Calculate y-intercepts
b <- mapply(function(x,y) calcB(-1/line[1], x, y), points$x[(start+1):(end-1)], points$y[(start+1):(end-1)])
#Calculate intercepts with with start-end line
ints <- sapply(b, function(x) intercept(line[1], line[2], -1/line[1], x), simplify=F)
#Calculate distances of points from line
distances <- mapply(function(x,y,z) distance(x[[1]], y, x[[2]], z), ints,
points$x[(start+1):(end-1)], points$y[(start+1):(end-1)])
#If any point greater than tolerance split at max distance point and apply to each side
if (any(distances >= tol)) {
keep[which.max(distances)+start] <<- T
#print(which.max(distances)+start)
simplify(start, which.max(distances)+start)
simplify(which.max(distances)+start, end)
}
}
}
#Start simplification
simplify(1, length(points$x))
#Return simplified points
list(x=points$x[keep], y=points$y[keep])
}
Try the shapefiles package in your browser
Any scripts or data that you put into this service are public.
shapefiles documentation built on Aug. 25, 2022, 5:05 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions dp change.id convert.to.simple convert.to.shapefile scaleXY add.xy write.shapefile read.shapefile write.dbf read.dbf write.shx calc.header read.shx write.shp read.shp
Documented in add.xy calc.header change.id convert.to.shapefile convert.to.simple dp read.dbf read.shapefile read.shp read.shx scaleXY write.dbf write.shapefile write.shp write.shx
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Shapefile Format - Read/Write shapefile format within R
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Created 2/13/03 Ben Stabler benjamin.stabler@odot.state.or.us
# Revised 6/11/03 Ben Stabler benjamin.stabler@odot.state.or.us
# Revised 7/7/03 Ben Stabler benjamin.stabler@odot.state.or.us
# Revised 7/23/03 Ben Stabler benjamin.stabler@odot.state.or.us
# Revised 8/10/05 Ben Stabler, benstabler@yahoo.com
# Reading of shape type 13 and 15 from Don MacQueen, macq@llnl.gov
# Conversion of simple polygons to shapefile format
# from Manuel Chirouze, Manuel.Chirouze@benfieldgroup.com
# Small bug fix to write.dbf to deal with 1 column white space correctly
# Revised 8/15/05 Ben Stabler, to use write.dbf and read.dbf from foreign library
# Added dp function as well
# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU General Public License
# as published by the Free Software Foundation; either version 2
# of the License, or (at your option) any later version.
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# DESCRIPTION
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#source("shapefiles.R")
#myshape <- read.shapefile("myshapewithnoextension")
#myshape$shp - the shp file
#myshape$shx - the index file
#myshape$dbf - the dbf file
#Each of the three lists of the myshape contain a header and some data. So there is...
#myshape$shp$header (the header info) and myshape$shp$shp (the geographic data)
#myshape$shx$header and myshape$shx$index (the index data)
#myshape$dbf$header and myshape$dbf$dbf (the dbf stored as a R data.frame)
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#Read in a point, line or polygon SHP file
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
read.shp <- function(shp.name) {
infile<-file(shp.name,"rb")
#Header
file.code <- readBin(infile,integer(),1, endian="big")
unused <- readBin(infile,integer(),5, endian="big")
file.length <- readBin(infile,integer(),1, endian="big")
file.version <- readBin(infile,integer(),1, endian="little")
shape.type <- readBin(infile,integer(),1, endian="little")
xmin <- readBin(infile,double(),1, endian="little")
ymin <- readBin(infile,double(),1, endian="little")
xmax <- readBin(infile,double(),1, endian="little")
ymax <- readBin(infile,double(),1, endian="little")
zmin <- readBin(infile,double(),1, endian="little")
zmax <- readBin(infile,double(),1, endian="little")
mmin <- readBin(infile,double(),1, endian="little")
mmax <- readBin(infile,double(),1, endian="little")
header <- list(file.code,file.length,file.version,shape.type,xmin,ymin,xmax,ymax,zmin,zmax,mmin,mmax)
names(header) <- c("file.code","file.length","file.version","shape.type","xmin","ymin","xmax","ymax","zmin","zmax","mmin","mmax")
rm(file.code,file.length,file.version,shape.type,xmin,ymin,xmax,ymax,zmin,zmax,mmin,mmax,unused)
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#Point Shapefile
if (header$shape.type==1) {
record <- integer()
shape.type <- integer()
x <- double()
y <- double()
record.num <- 1
while (length(record.num)!=0) {
#Record Number
record.num <- readBin(infile,integer(),1, 4, endian="big")
if (length(record.num)==0) break
record <- c(record,record.num)
#Content Length
content.length <- readBin(infile,integer(),1, 4, endian="big")
#Shape Type (should be 1)
temp.type <- readBin(infile,integer(),1, 4, endian="little")
shape.type <- c(shape.type, temp.type)
if (temp.type==1) {
#X and Y values
x <- c(x,readBin(infile,double(),1, 8, endian="little"))
y <- c(y,readBin(infile,double(),1, 8, endian="little"))
#To read in extra data created by GeoMedia shapefile export
if (content.length>10) {
temp <- readBin(infile,integer(), ((content.length-10)/2), 4, endian="big")
}
}
if (temp.type==0) {
#Null shape type has no data
x <- c(x,NA)
y <- c(y,NA)
}
}
point.data <- cbind(record=record, x=x, y=y, shape.type=shape.type)
close(infile)
}
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#PolyLine, Polygon, PolyLineZ, and PolygonZ Shapefile
if (header$shape.type==3 || header$shape.type==5 || header$shape.type == 13 || header$shape.type == 15) {
shape <- list()
record.num <- 1
while (length(record.num)!=0) {
#Record Number
record.num <- readBin(infile,integer(),1, endian="big")
if (length(record.num)==0) break
#Content Length
content.length <- readBin(infile,integer(),1, endian="big")
#Shape Type
shape.type <- readBin(infile, integer(), 1, endian="little")
if (shape.type==3 || shape.type==5) {
box <- readBin(infile, double(), 4, endian="little")
names(box) <- c("xmin","ymin","xmax","ymax")
num.parts <- readBin(infile, integer(), 1, endian="little")
num.points <- readBin(infile, integer(), 1, endian="little")
parts <- readBin(infile, integer(), num.parts, endian="little")
points <- readBin(infile, double(), num.points*2, endian="little")
X <- points[seq(1,(num.points*2),by=2)]
Y <- points[seq(2,(num.points*2),by=2)]
points <- data.frame(X,Y)
}
#Type 13 polylineZ (and 15 polygonZ since it is the same as 13)
if (shape.type==13 || shape.type==15) {
box <- readBin(infile, double(), 4, endian = "little")
names(box) <- c("xmin", "ymin", "xmax", "ymax")
num.parts <- readBin(infile, integer(), 1, endian = "little")
num.points <- readBin(infile, integer(), 1, endian = "little")
parts <- readBin(infile, integer(), num.parts, endian = "little")
points <- readBin(infile, double(), num.points * 2, endian = "little")
X <- points[seq(1, (num.points * 2), by = 2)]
Y <- points[seq(2, (num.points * 2), by = 2)]
#Additional information for type 13 and 15
zrange <- readBin(infile,double(),2,endian='little')
Z <- readBin(infile,double(),num.points,endian='little')
mrange <- readBin(infile,double(),2,endian='little')
M <- readBin(infile,double(),num.points,endian='little')
points <- data.frame(X, Y, Z, M)
}
if (shape.type==0) {
#Null shape type has no data
box <- rep(NA,4)
num.parts <- NA
num.points <- NA
parts <- NA
points <- data.frame(X=NA,Y=NA)
}
shape.info <- list(record=record.num, content.length=content.length, shape.type=shape.type, box=box, num.parts=num.parts, num.points=num.points, parts=parts, points=points)
shape <- c(shape,list(shape.info))
}
close(infile)
}
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if (header$shape.type==1) {
return.val <- point.data
}
if (header$shape.type==3 || header$shape.type==5 || header$shape.type == 13 || header$shape.type == 15) {
return.val <- shape
}
list(shp=return.val,header=header)
}
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#Write out a point, line or polygon SHP file
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
write.shp <- function(shp, out.name) {
outfile<-file(out.name,"wb")
#Header
writeBin(as.integer(9994), outfile, endian="big")
writeBin(as.integer(rep(0,5)), outfile, endian="big")
writeBin(as.integer(shp$header$file.length), outfile, endian="big")
writeBin(as.integer(1000), outfile, endian="little")
writeBin(as.integer(shp$header$shape.type), outfile, endian="little")
writeBin(as.double(shp$header$xmin), outfile, endian="little")
writeBin(as.double(shp$header$ymin), outfile, endian="little")
writeBin(as.double(shp$header$xmax), outfile, endian="little")
writeBin(as.double(shp$header$ymax), outfile, endian="little")
writeBin(as.double(shp$header$zmin), outfile, endian="little")
writeBin(as.double(shp$header$zmax), outfile, endian="little")
writeBin(as.double(shp$header$mmin), outfile, endian="little")
writeBin(as.double(shp$header$mmax), outfile, endian="little")
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if (shp$header$shape.type==1) {
for (record in 1:nrow(shp$shp)) {
if (shp$shp[record,4]==1) {
#Record Number
writeBin(as.integer(shp$shp[record,1]), outfile, endian="big")
#Content Length
writeBin(as.integer(10), outfile, endian="big")
#Shape Type
writeBin(as.integer(1), outfile, endian="little")
#X - For Points only
writeBin(as.double(shp$shp[record,2]), outfile, endian="little")
#Y - For Points only
writeBin(as.double(shp$shp[record,3]), outfile, endian="little")
}
if (shp$shp[record,4]==0) {
#Record Number
writeBin(as.integer(shp$shp[record,1]), outfile, endian="big")
#Content Length
writeBin(as.integer(4), outfile, endian="big")
#Shape Type
writeBin(as.integer(0), outfile, endian="little")
}
}
}
if (shp$header$shape.type==3 || shp$header$shape.type==5 || shp$header$shape.type == 13 || shp$header$shape.type == 15) {
for (record in 1:length(shp$shp)) {
if (shp$shp[[record]]$shape.type==3 || shp$shp[[record]]$shape.type==5) {
#Record Number
writeBin(as.integer(shp$shp[[record]]$record), outfile, endian="big")
#Content Length
writeBin(as.integer(shp$shp[[record]]$content.length), outfile, endian="big")
#Shape Type
writeBin(as.integer(shp$shp[[record]]$shape.type), outfile, endian="little")
#Bounding Box
writeBin(as.double(shp$shp[[record]]$box[1]), outfile, endian="little")
writeBin(as.double(shp$shp[[record]]$box[2]), outfile, endian="little")
writeBin(as.double(shp$shp[[record]]$box[3]), outfile, endian="little")
writeBin(as.double(shp$shp[[record]]$box[4]), outfile, endian="little")
#Number of parts (segments)
writeBin(as.integer(shp$shp[[record]]$num.parts), outfile, endian="little")
#Number of total points
writeBin(as.integer(shp$shp[[record]]$num.points), outfile, endian="little")
#Parts
writeBin(as.integer(shp$shp[[record]]$parts), outfile, endian="little")
#Points - merge the X and Y points into one vector
point.stream <- rep(0,length(shp$shp[[record]]$points$X)*2)
point.stream[seq(1,length(shp$shp[[record]]$points$X)*2,by=2)] <- shp$shp[[record]]$points$X
point.stream[seq(2,length(shp$shp[[record]]$points$Y)*2,by=2)] <- shp$shp[[record]]$points$Y
writeBin(as.double(point.stream), outfile, endian="little")
#Need to check to make sure first and last vertex are the same for Polygon
}
#Type 13 polylineZ (and 15 polygonZ since it is the same as 13)
if (shp$shp[[record]]$shape.type==13 || shp$shp[[record]]$shape.type==15) {
#Record Number
writeBin(as.integer(shp$shp[[record]]$record), outfile, endian="big")
#Content Length
writeBin(as.integer(shp$shp[[record]]$content.length), outfile, endian="big")
#Shape Type
writeBin(as.integer(shp$shp[[record]]$shape.type), outfile, endian="little")
#Bounding Box
writeBin(as.double(shp$shp[[record]]$box[1]), outfile, endian="little")
writeBin(as.double(shp$shp[[record]]$box[2]), outfile, endian="little")
writeBin(as.double(shp$shp[[record]]$box[3]), outfile, endian="little")
writeBin(as.double(shp$shp[[record]]$box[4]), outfile, endian="little")
#Number of parts (segments)
writeBin(as.integer(shp$shp[[record]]$num.parts), outfile, endian="little")
#Number of total points
writeBin(as.integer(shp$shp[[record]]$num.points), outfile, endian="little")
#Parts
writeBin(as.integer(shp$shp[[record]]$parts), outfile, endian="little")
#Points - merge the X and Y points into one vector
point.stream <- rep(0,length(shp$shp[[record]]$points$X)*2)
point.stream[seq(1,length(shp$shp[[record]]$points$X)*2,by=2)] <- shp$shp[[record]]$points$X
point.stream[seq(2,length(shp$shp[[record]]$points$Y)*2,by=2)] <- shp$shp[[record]]$points$Y
writeBin(as.double(point.stream), outfile, endian="little")
#Need to check to make sure first and last vertex are the same for Polygon
#Write out Z Range and Z values
minZ <- min(shp$shp[[record]]$points$Z)
maxZ <- max(shp$shp[[record]]$points$Z)
writeBin(as.double(minZ), outfile, endian='little')
writeBin(as.double(maxZ), outfile, endian='little')
writeBin(as.double(shp$shp[[record]]$points$Z), outfile, endian='little')
#Write out M Range and M values
minM <- min(shp$shp[[record]]$points$M)
maxM <- max(shp$shp[[record]]$points$M)
writeBin(as.double(minM), outfile, endian='little')
writeBin(as.double(maxM), outfile, endian='little')
writeBin(as.double(shp$shp[[record]]$points$M), outfile, endian='little')
}
if (shp$shp[[record]]$shape.type==0) {
#Record Number
writeBin(as.integer(shp$shp[[record]]$record), outfile, endian="big")
#Content Length
writeBin(as.integer(shp$shp[[record]]$content.length), outfile, endian="big")
#Shape Type
writeBin(as.integer(shp$shp[[record]]$shape.type), outfile, endian="little")
}
}
}
close(outfile)
}
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#Read in the SHX file (index)
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
read.shx <- function(shx.name) {
infile<-file(shx.name,"rb")
#Header
file.code <- readBin(infile,integer(),1, endian="big")
unused <- readBin(infile,integer(),5, endian="big")
file.length <- readBin(infile,integer(),1, endian="big")
file.version <- readBin(infile,integer(),1, endian="little")
shape.type <- readBin(infile,integer(),1, endian="little")
xmin <- readBin(infile,double(),1, endian="little")
ymin <- readBin(infile,double(),1, endian="little")
xmax <- readBin(infile,double(),1, endian="little")
ymax <- readBin(infile,double(),1, endian="little")
zmin <- readBin(infile,double(),1, endian="little")
zmax <- readBin(infile,double(),1, endian="little")
mmin <- readBin(infile,double(),1, endian="little")
mmax <- readBin(infile,double(),1, endian="little")
header <- list(file.code,file.length,file.version,shape.type,xmin,ymin,xmax,ymax,zmin,zmax,mmin,mmax)
names(header) <- c("file.code","file.length","file.version","shape.type","xmin","ymin","xmax","ymax","zmin","zmax","mmin","mmax")
records <- (file.length-50)/4
#Record Data
index.data <- matrix(0,records,2)
colnames(index.data) <- c("Offset","Length")
for (record in 1:records) {
#First record offset is 50
index.data[record,1] <- readBin(infile,integer(),1, endian="big")
#Content length
index.data[record,2] <- readBin(infile,integer(),1, endian="big")
}
close(infile)
list(index=index.data,header=header)
}
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Function to calculate header information for a point file
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Used to clean up GeoMedia shapefile export files
calc.header <- function(shapefile) {
#1) Calculate shp header info
#Calculate number of bytes for null record headers
num.null.record.bytes <- length(shapefile$shp$shp[,2][is.na(shapefile$shp$shp[,2])]) * 4
#Calculate number of bytes for complete record record headers
num.record.header.bytes <- length(shapefile$shp$shp[,2]) * 8
#Calculate number of bytes for records (both geoMedia and generic shapefile)
#Should be value 20 (56 for GeoMedia files)
num.record.bytes <- (length(shapefile$shp$shp[,2]) - length(shapefile$shp$shp[,2][is.na(shapefile$shp$shp[,2])])) * 20
#Sum the bytes to get total bytes (100 = header bytes)
total.bytes <- 100 + sum(num.null.record.bytes, num.record.header.bytes, num.record.bytes)
#Divide by 2 to get number of 16-bit words
file.length <- total.bytes / 2
#Replace file.length in shapefile header with new file.length calculation
shapefile$shp$header$file.length <- file.length
#Replace lengths of 28 with 10 for shx from GeoMedia
#2) Calculate shx info
if (shapefile$shx$index[1,2]==28) {
shapefile$shx$index[,2][shapefile$shx$index[,2]==28] <- 10
}
#Set the content.length to 14 in order to calculate offset
shapefile$shx$index[,2][shapefile$shx$index[,2]==10]<-14
#Set null value content.lengths from 2 to 6 to calculate offsets
shapefile$shx$index[,2][shapefile$shx$index[,2]==2]<-6
#Calculate new offset values for each record from beginning of shp file
new.offset <- cumsum(c(50,shapefile$shx$index[,2]))
#Replace the existing offsets with the new ones
shapefile$shx$index[,1] <- new.offset[1:(length(new.offset)-1)]
#Set the content.lenghts back to their respective values
shapefile$shx$index[,2][shapefile$shx$index[,2]==14]<-10
shapefile$shx$index[,2][shapefile$shx$index[,2]==6]<-2
shapefile
}
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#Write out SHX file (index)
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
write.shx <- function(shx, out.name) {
outfile<-file(out.name,"wb")
#Header
writeBin(as.integer(9994), outfile, endian="big")
writeBin(as.integer(rep(0,5)), outfile, endian="big")
writeBin(as.integer(shx$header$file.length), outfile, endian="big")
writeBin(as.integer(1000), outfile, endian="little")
writeBin(as.integer(shx$header$shape.type), outfile, endian="little")
writeBin(as.double(shx$header$xmin), outfile, endian="little")
writeBin(as.double(shx$header$ymin), outfile, endian="little")
writeBin(as.double(shx$header$xmax), outfile, endian="little")
writeBin(as.double(shx$header$ymax), outfile, endian="little")
writeBin(as.double(shx$header$zmin), outfile, endian="little")
writeBin(as.double(shx$header$zmax), outfile, endian="little")
writeBin(as.double(shx$header$mmin), outfile, endian="little")
writeBin(as.double(shx$header$mmax), outfile, endian="little")
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
for (record in 1:nrow(shx$index)) {
#Record Offset
writeBin(as.integer(shx$index[record,1]), outfile, endian="big")
#Content Length
writeBin(as.integer(shx$index[record,2]), outfile, endian="big")
}
close(outfile)
}
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#Read DBF format
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
read.dbf <- function(dbf.name, header=FALSE) {
#Load foreign package for read.dbf
## library(foreign)
#Read header if needed
if(header) {
infile<-file(dbf.name,"rb")
#Header
file.version <- readBin(infile,integer(), 1, size=1, endian="little")
file.year <- readBin(infile,integer(), 1, size=1, endian="little")
file.month <- readBin(infile,integer(), 1, size=1, endian="little")
file.day <- readBin(infile,integer(), 1, size=1, endian="little")
num.records <- readBin(infile,integer(), 1, size=4, endian="little")
header.length <- readBin(infile,integer(), 1, size=2, endian="little")
record.length <- readBin(infile,integer(), 1, size=2, endian="little")
file.temp <- readBin(infile,integer(), 20, size=1, endian="little")
header <- list(file.version,file.year, file.month, file.day, num.records, header.length, record.length)
names(header) <- c("file.version","file.year","file.month","file.day","num.records","header.length","record.length")
rm(file.version,file.year, file.month, file.day, num.records, header.length, record.length)
#Calculate the number of fields
num.fields <- (header$header.length-32-1)/32
field.name <- NULL
field.type <- NULL
field.length <- NULL
field.decimal <- NULL
#Field Descriptions (32 bytes each)
for (i in 1:num.fields) {
field.name.test <- readBin(infile,character(), 1, size=10, endian="little")
field.name <- c(field.name,field.name.test)
if (nchar(field.name.test)!=10) {
file.temp <- readBin(infile,integer(), 10-(nchar(field.name.test)), 1, endian="little")
}
field.type <- c(field.type,readChar(infile, 1))
file.temp <- readBin(infile,integer(), 4, 1, endian="little")
field.length <- c(field.length,readBin(infile,integer(), 1, 1, endian="little"))
field.decimal <- c(field.decimal, readBin(infile,integer(), 1, 1, endian="little"))
file.temp <- readBin(infile,integer(), 14, 1, endian="little")
}
#Create a table of the field info
fields <- data.frame(NAME=field.name,TYPE=field.type,LENGTH=field.length,DECIMAL=field.decimal)
#Set all fields with length<0 equal to correct number of characters
fields$LENGTH[fields$LENGTH<0]<-(256+fields$LENGTH[fields$LENGTH<0])
#Read in end of attribute descriptions terminator - should be integer value 13
file.temp <- readBin(infile,integer(), 1, 1, endian="little")
#Increase the length of field 1 by one to account for the space at the beginning of each record
fields$LENGTH[1]<-fields$LENGTH[1]+1
#Add fields to the header list
header <- c(header,fields=NULL)
header$fields <- fields
#Close file connection
close(infile)
}
#Read in dbf
dbf <- foreign::read.dbf(dbf.name)
#Return the dbf as a list with a data.frame and a header list
list(dbf=dbf, header=header)
}
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#Write out DBF format
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
write.dbf <- function(dbf, out.name, arcgis=FALSE) {
#Load foreign package for write.dbf
## library(foreign)
#If output intended for ArcGIS, replace "." with "_" for column names in header
if (arcgis==T) {
colnames(dbf$dbf) <- gsub("\\.","_",colnames(dbf$dbf))
}
#Write out dbf
foreign::write.dbf(dbf$dbf, out.name)
}
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#Function to Read in a Shapefile
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
read.shapefile <- function(shape.name) {
shp.data <- read.shp(paste(shape.name, ".shp", sep=""))
shx.data <- read.shx(paste(shape.name, ".shx", sep=""))
dbf.data <- get("read.dbf","package:shapefiles")(paste(shape.name, ".dbf", sep=""))
shapefile <- list(shp=shp.data,shx=shx.data,dbf=dbf.data)
shapefile
}
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#Function to Write out a Shapefile
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
write.shapefile <- function(shapefile, out.name, arcgis=FALSE) {
write.shp(shapefile$shp, paste(out.name, ".shp", sep=""))
write.shx(shapefile$shx, paste(out.name, ".shx", sep=""))
get("write.dbf","package:shapefiles")(shapefile$dbf, paste(out.name, ".dbf", sep=""), arcgis)
}
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#Function to Add the X and Y Coordinates to the DBF of a shapefile
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
add.xy <- function(shapefile) {
if (shapefile$shp$header$shape.type !=1) {
stop ("Must be point shapefile")
}
shapefile$dbf[[1]]$XCOORD <- shapefile$shp[[1]][,2]
shapefile$dbf[[1]]$YCOORD <- shapefile$shp[[1]][,3]
shapefile
}
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#Function to scale the X and Y Coordinates of the shapefile
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
scaleXY <- function(shapefile, scale.factor) {
#Scale the bounding box of the shapefile by the scale factor
shapefile$shp$header$xmin <- shapefile$shp$header$xmin/scale.factor
shapefile$shp$header$ymin <- shapefile$shp$header$ymin/scale.factor
shapefile$shp$header$xmax <- shapefile$shp$header$xmax/scale.factor
shapefile$shp$header$ymax <- shapefile$shp$header$ymax/scale.factor
#Scale the bounding box of the shx file by the scale factor
shapefile$shx$header$xmin <- shapefile$shx$header$xmin/scale.factor
shapefile$shx$header$ymin <- shapefile$shx$header$ymin/scale.factor
shapefile$shx$header$xmax <- shapefile$shx$header$xmax/scale.factor
shapefile$shx$header$ymax <- shapefile$shx$header$ymax/scale.factor
#Scale the X and Y if a point shape
if (shapefile$shp$header$shape.type == 1) {
shapefile$shp[[1]][,2] <- shapefile$shp[[1]][,2]/scale.factor
shapefile$shp[[1]][,3] <- shapefile$shp[[1]][,3]/scale.factor
}
#Scale the X and Y point values and bounding box if a line or polyon shape
if (shapefile$shp$header$shape.type==3 || shapefile$shp$header$shape.type==5) {
for (shape in 1:length(shapefile$shp$shp)) {
shapefile$shp[[1]][[shape]]$points[,1] <- shapefile$shp[[1]][[shape]]$points[,1]/scale.factor
shapefile$shp[[1]][[shape]]$points[,2] <- shapefile$shp[[1]][[shape]]$points[,2]/scale.factor
shapefile$shp[[1]][[shape]]$box <- shapefile$shp[[1]][[shape]]$box/scale.factor
}
}
shapefile
}
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#Converts a simple point, polyLine, or polygon data frame into a shapefile.
#The shpTable data frame must have three columns in this order: Id, X, and Y.
#The attTable data frame must have a field to join on - identified by field
#The type argument is either 1 (point), 3 (polyline) or 5 (polygon).
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
convert.to.shapefile <- function(shpTable, attTable, field, type) {
#Check to ensure field in attTable
if(!(field %in% colnames(attTable))) {
stop("Field not in attTable column names")
}
#Check to ensure the same Id set in shpTable and attTable
if(length(unique(shpTable[,1])) != length(attTable[,field])) {
stop("Different number of unique Ids in shpTable versus attTable")
}
if(any(sort(unique(shpTable[,1])) != sort(attTable[,field]))) {
stop("Id set in shpTable not the same as Id set in attTable")
}
#Build dbf - no longer uses header since "foreign.write.dbf" does not
dbf <- list(dbf=attTable,header=NULL)
#Build shp file
#shp list
shpTable[,c(2,3)] <- round(shpTable[,c(2,3)], 5)
colnames(shpTable) <- c("Id","X","Y")
v.content.length <- vector()
#point
if(type == 1) {
#loop for each entry in the dbf and build a shape
for (i in 1:nrow(attTable)) {
#Get shape entries for attTable Id i
shpTable.i <- shpTable[which(shpTable[,1] == attTable[i,field]), -1]
p.content.length <- 10 #shape.type and 2 doubles
v.content.length[i] <- p.content.length
if(i == 1) {
shp <- data.frame(
record = i,
x = shpTable.i[,1],
y = shpTable.i[,2],
shape.type = type)
} else {
shp2 <- data.frame(
record = i,
x = shpTable.i[,1],
y = shpTable.i[,2],
shape.type = type)
shp <- rbind(shp, shp2)
rm(shp2)
}
}
}
#polyLine and polygon shp
if(type == 3 || type == 5) {
#shp file
shp<-list()
#loop for each entry in the dbf and build a shape
for (i in 1:nrow(attTable)) {
#Get shape entries for attTable Id i
shpTable.i <- shpTable[which(shpTable[,1] == attTable[i,field]), -1]
num.points <- nrow(shpTable.i)
v.box <- c(min(shpTable.i[,1]),min(shpTable.i[,2]),max(shpTable.i[,1]),max(shpTable.i[,2]))
names(v.box) <- c("xmin","ymin","xmax","ymax")
p.content.length <- 24 + (num.points * 8)
v.content.length[i] <- p.content.length
shp[[i]]<-list(record = i,
content.length = p.content.length,
shape.type = type,
box = v.box,
num.parts = 1, #Only 1 part
num.points = nrow(shpTable.i),
parts = 0,
points = shpTable.i)
}
}
#Collect up all shp info
header<-list(file.code = 9994,
file.length = 50 + sum(v.content.length + 4),
file.version = 1000,
shape.type = type,
xmin = min(shpTable[,2]),
ymin = min(shpTable[,3]),
xmax = max(shpTable[,2]),
ymax = max(shpTable[,3]),
zmin = 0,
zmax = 0,
mmin = 0,
mmax = 0)
shp <- list(shp=shp, header=header)
#Build shx file
#shx list
shx<-list(index = cbind(Offset=c(50 ,(cumsum(v.content.length + 4) + 50)[-length(v.content.length)]),
Length = v.content.length),
header = shp$header)
shx$header$file.length <- 50 + 4 * nrow(shx$index)
#Return shapefile
shp.file <- list(shp=shp, shx=shx, dbf=dbf)
shp.file
}
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#Convert read.shp result to simple data.frame format (Id, X, Y)
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
convert.to.simple <- function(shp) {
#Determine type of shp
if(shp$header$shape.type == 1) {
#Get points
allPoints <- shp$shp[,c("record","x","y")]
colnames(allPoints) <- c("Id", "X", "Y")
} else {
#Get points
allPoints <- do.call(rbind, lapply(shp$shp, function(x) x$points[,c("X","Y")]))
#Calculate index and add to points
indexes <- 1:length(shp$shp)
repNum <- lapply(shp$shp, function(x) nrow(x$points))
allPoints <- cbind(Id=rep(indexes, repNum), allPoints)
}
#Return allPoints
allPoints
}
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#Change Index to a Field from the DBF
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#Changes the Index field of the shpTable to the new vector of fields
#Will replace the records in order, so Index 1 = element 1 in the new vector
change.id <- function(shpTable, newFieldAsVector) {
#Check to ensure the same Id set in shpTable and attTable
if(length(unique(shpTable[,1])) != length(newFieldAsVector)) {
stop("Different number of unique Ids in shpTable versus new field")
}
#Create new vector
count <- table(shpTable[,1])
shpTable[,1] <- rep(newFieldAsVector, count)
#Return result
shpTable
}
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#Douglas-Peucker Polyline Simplification
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#Ben Stabler, bstabler@ptvamerica.com, March 2005
#Douglas, D. and Peucker, T. (1973). "Algorithms for
#the reduction of the number of points required to
#represent a digitized line or its caricature."
#The Canadian Cartographer 10(2). 112-122.
#Currently uses the line, not the line segment to
#determine the distance of the points from the line.
#See http://www.lgc.com/resources/Doug_Peucker.pdf
#for more information. This can result in the
#omission of extreme "outlier-like" points. Try:
#points<-list(x=c(200,0,400),y=c(0,20,0))
#plot(points, type="l")
#lines(dp(points, 50), type="l", col="blue")
#Simple Example
#x <- c(5,3,4,1,8,9,10,11)
#y <- c(6,4,2,1,1,5,2,3)
#points <- list(x=x,y=y)
#plot(points, type="l")
#lines(dp(points, 2), type="l", col="blue")
####################################################
dp <- function(points, tolerance) {
#Convert to lowercase
names(points) <- tolower(names(points))
#Calculate distance between two points
distance <- function(x1,x2,y1,y2) {
sqrt((x2-x1)^2 + (y2-y1)^2)
}
#Calculate equation of a line from two points
equationOfLine <- function(x1,x2,y1,y2) {
slope <- (y2-y1)/(x2-x1)
b <- y1 - slope * x1
c(slope,b)
}
#Calculate y-intercept from a point and a slope
calcB <- function(slope, px, py) {
py - (slope * px)
}
#Calculate intercept of two lines from their slope and y-intercept
intercept <- function(s1,b1,s2,b2) {
x1 <- (b2-b1)/(s1-s2)
y1 <- s1 * x1 + b1
c(x1,y1)
}
#Setup vector to mark points to keep
keep <- rep(F, length(points$x))
keep[1] <- T
keep[length(keep)] <- T
#Function definition to simplify points
simplify <- function(start, end, tol=tolerance) {
#Calculate intermediate point distances
if (length(points$x[start:end]) > 2) {
#Avoid Inf slope
if( points$x[start] == points$x[end] ) { points$x[start] <- points$x[start] - 0.0000001 }
if( points$y[start] == points$y[end]) { points$y[start] <- points$y[start] - 0.0000001 }
#Calculate equation of line of middle points
line <- equationOfLine( points$x[start], points$x[end], points$y[start], points$y[end])
#Calculate y-intercepts
b <- mapply(function(x,y) calcB(-1/line[1], x, y), points$x[(start+1):(end-1)], points$y[(start+1):(end-1)])
#Calculate intercepts with with start-end line
ints <- sapply(b, function(x) intercept(line[1], line[2], -1/line[1], x), simplify=F)
#Calculate distances of points from line
distances <- mapply(function(x,y,z) distance(x[[1]], y, x[[2]], z), ints,
points$x[(start+1):(end-1)], points$y[(start+1):(end-1)])
#If any point greater than tolerance split at max distance point and apply to each side
if (any(distances >= tol)) {
keep[which.max(distances)+start] <<- T
#print(which.max(distances)+start)
simplify(start, which.max(distances)+start)
simplify(which.max(distances)+start, end)
}
}
}
#Start simplification
simplify(1, length(points$x))
#Return simplified points
list(x=points$x[keep], y=points$y[keep])
}
Try the shapefiles package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.