Tools for reading and handling spatial objects

Tools for reading and handling spatial objects [R logo]

Help pages for package `maptools' version 0.6-9

A C D G H I K L M N O P R S T U W

-- A --

ArcObj2SLDF | Converter functions to build SpatialLinesDataFrame objects
as.im.SpatialGridDataFrame | coercion between sp objects and spatstat objects
as.owin.SpatialGridDataFrame | coercion between sp objects and spatstat objects
as.owin.SpatialPixelsDataFrame | coercion between sp objects and spatstat objects
as.owin.SpatialPolygons | coercion between sp objects and spatstat objects
as.ppp.SpatialGridDataFrame | coercion between sp objects and spatstat objects
as.ppp.SpatialPoints | coercion between sp objects and spatstat objects
as.ppp.SpatialPointsDataFrame | coercion between sp objects and spatstat objects
as.SpatialGridDataFrame.im | coercion between sp objects and spatstat objects
as.SpatialGridDataFrame.ppp | coercion between sp objects and spatstat objects
as.SpatialPoints.ppp | coercion between sp objects and spatstat objects
as.SpatialPointsDataFrame.ppp | coercion between sp objects and spatstat objects

-- C --

checkPolygonsHoles | Check holes in Polygons objects
coerce,SpatialGridDataFrame,im-method | coercion between sp objects and spatstat objects
coerce,SpatialGridDataFrame,owin-method | coercion between sp objects and spatstat objects
coerce,SpatialPixelsDataFrame,owin-method | coercion between sp objects and spatstat objects
coerce,SpatialPoints,ppp-method | coercion between sp objects and spatstat objects
coerce,SpatialPointsDataFrame,ppp-method | coercion between sp objects and spatstat objects
coerce,SpatialPolygons,owin-method | coercion between sp objects and spatstat objects
ContourLines2SLDF | Converter functions to build SpatialLinesDataFrame objects
convert.pl | Create polygon lists and bounding boxes from imported shapefiles
crepuscule | Methods for sun ephemerides calculations
crepuscule,matrix,POSIXct-method | Methods for sun ephemerides calculations
crepuscule,SpatialPoints,POSIXct-method | Methods for sun ephemerides calculations
crepuscule-methods | Methods for sun ephemerides calculations

-- D --

dotsInPolys Put dots in polygons

-- G --

get.Pcent | Polygon centroids
getinfo.shape | Read shapefile into Map object
GE_SpatialGrid | Create SpatialGrid for PNG output to GE
gpcholes | Hisaji Ono's lake/hole problem
gzAzimuth | Find azimuth for geographical coordinates

-- H --

hole1pl | Hisaji Ono's lake/hole problem
hole2pl | Hisaji Ono's lake/hole problem

-- I --

im-class Virtual class "ppp"

-- K --

kmlOverlay Create and write KML file for PNG image overlay

-- L --

leglabs Plot polygons

-- M --

Map2bbs | Create polygon lists and bounding boxes from imported shapefiles
Map2lines | Create polygon lists and bounding boxes from imported shapefiles
Map2points | Create polygon lists and bounding boxes from imported shapefiles
Map2poly | Create polygon lists and bounding boxes from imported shapefiles
Map2poly1 | Create polygon lists and bounding boxes from imported shapefiles
map2SpatialLines | Convert map objects to sp classes
map2SpatialPolygons | Convert map objects to sp classes
MapGen2SL | Converter functions to build SpatialLinesDataFrame objects
MapShapeIds | Create polygon lists and bounding boxes from imported shapefiles
maptools | Report version information and changes

-- N --

nowrapRecenter | Break polygons at meridian for recentering
nowrapSpatialPolygons | Break polygons at meridian for recentering

-- O --

owin-class Virtual class "ppp"

-- P --

pal2SpatialPolygons | Making SpatialPolygons objects from RArcInfo input
plot.Map | Plot a Map object (deprecated)
plot.polylist | Plot polygons
plotpolys | Plot polygons
ppp-class | Virtual class "ppp"
print.shapehead | Read shapefile into Map object
pruneMap | Convert map objects to sp classes

-- R --

read.shape | Read shapefile into Map object
readAsciiGrid | read/write to/from (ESRI) asciigrid format
readGPS | GPSbabel read interface
readShapeLines | Read arc shape files into SpatialLinesDataFrame objects
readShapePoints | Read points shape files into SpatialPointsDataFrame objects
readShapePoly | Read polygon shape files into SpatialPolygonsDataFrame objects
Rgshhs | Read GSHHS polygons into SpatialPolygons object

-- S --

shape2bbs | Create polygon lists and bounding boxes from imported shapefiles
shape2lines | Create polygon lists and bounding boxes from imported shapefiles
shape2points | Create polygon lists and bounding boxes from imported shapefiles
shape2poly | Create polygon lists and bounding boxes from imported shapefiles
solarnoon | Methods for sun ephemerides calculations
solarnoon,matrix,POSIXct-method | Methods for sun ephemerides calculations
solarnoon,SpatialPoints,POSIXct-method | Methods for sun ephemerides calculations
solarnoon-methods | Methods for sun ephemerides calculations
solarpos | Methods for sun ephemerides calculations
solarpos,matrix,POSIXct-method | Methods for sun ephemerides calculations
solarpos,SpatialPoints,POSIXct-method | Methods for sun ephemerides calculations
solarpos-methods | Methods for sun ephemerides calculations
sp2Mondrian | write map data for Mondrian
sp2tmap | Convert SpatialPolygons object for Stata tmap command
sp2WB | Export SpatialPolygons object as S-Plus map for WinBUGS
SpatialLines2PolySet | Convert sp line and polygon objects to PBSmapping PolySet objects
SpatialPolygons2PolySet | Convert sp line and polygon objects to PBSmapping PolySet objects
spCbind | cbind for spatial objects
spCbind,SpatialLinesDataFrame,data.frame-method | cbind for spatial objects
spCbind,SpatialLinesDataFrame,vector-method | cbind for spatial objects
spCbind,SpatialPointsDataFrame,data.frame-method | cbind for spatial objects
spCbind,SpatialPointsDataFrame,vector-method | cbind for spatial objects
spCbind,SpatialPolygonsDataFrame,data.frame-method | cbind for spatial objects
spCbind,SpatialPolygonsDataFrame,vector-method | cbind for spatial objects
spCbind-methods | cbind for spatial objects
spChFIDs | change feature IDs in spatial objects
spChFIDs,SpatialLines,character-method | change feature IDs in spatial objects
spChFIDs,SpatialLinesDataFrame,character-method | change feature IDs in spatial objects
spChFIDs,SpatialPolygons,character-method | change feature IDs in spatial objects
spChFIDs,SpatialPolygonsDataFrame,character-method | change feature IDs in spatial objects
spChFIDs-methods | change feature IDs in spatial objects
spRbind | rbind for spatial objects
spRbind,SpatialLines,SpatialLines-method | rbind for spatial objects
spRbind,SpatialLinesDataFrame,SpatialLinesDataFrame-method | rbind for spatial objects
spRbind,SpatialPoints,SpatialPoints-method | rbind for spatial objects
spRbind,SpatialPointsDataFrame,SpatialPointsDataFrame-method | rbind for spatial objects
spRbind,SpatialPolygons,SpatialPolygons-method | rbind for spatial objects
spRbind,SpatialPolygonsDataFrame,SpatialPolygonsDataFrame-method | rbind for spatial objects
spRbind-methods | rbind for spatial objects
subset.polylist | Subset polygon list objects
sunriset | Methods for sun ephemerides calculations
sunriset,matrix,POSIXct-method | Methods for sun ephemerides calculations
sunriset,SpatialPoints,POSIXct-method | Methods for sun ephemerides calculations
sunriset-methods | Methods for sun ephemerides calculations
symbolsInPolys | Place grids of points over polygons

-- T --

trackAzimuth Find azimuth for geographical coordinates

-- U --

unionSpatialPolygons Aggregate Polygons in a SpatialPolygons object

-- W --

write.linelistShape | Write a arc-type shapefile
write.pointShape | Write a point-type shapefile
write.polylistShape | Write a polygon-type shapefile
writeAsciiGrid | read/write to/from (ESRI) asciigrid format
writeLinesShape | Read arc shape files into SpatialLinesDataFrame objects
writePointsShape | Read points shape files into SpatialPointsDataFrame objects
writePolyShape | Read polygon shape files into SpatialPolygonsDataFrame objects


coercion between sp objects and spatstat objects

Description

S4-style as() coercion works from objects of S4 sp classes to spatstat S3 classes, but not the reverse. Coercion to sp classes uses direct function calls.

Usage

as.SpatialPoints.ppp(from)
as.SpatialPointsDataFrame.ppp(from)
as.SpatialGridDataFrame.ppp(from)
as.SpatialGridDataFrame.im(from)

Arguments

from spatstat object of class ppp or im

Methods

coerce
signature(from = "SpatialPoints", to = "ppp")
coerce
signature(from = "SpatialPointsDataFrame", to = "ppp")
coerce
signature(from = "SpatialGridDataFrame", to = "ppp")
coerce
signature(from = "SpatialPolygons", to = "owin")
coerce
signature(from = "SpatialPixelsDataFrame", to = "owin")
coerce
signature(from = "SpatialGridDataFrame", to = "owin")
coerce
signature(from = "SpatialGridDataFrame", to = "im")

Author(s)

Edzer Pebesma e.pebesma@geo.uu.nl, Roger Bivand

Examples

library(spatstat)
data(meuse)
coordinates(meuse) = ~x+y
zn1 <- as(meuse["zinc"], "ppp")
zn1
plot(zn1)
as(as(meuse, "SpatialPoints"), "ppp")
data(meuse.grid)
gridded(meuse.grid) = ~x+y
mg_owin <- as(meuse.grid, "owin")
zn1a <- ppp(x=zn1$x, y=zn1$y, marks=zn1$marks, window=mg_owin)
zn1a
plot(zn1a)
rev_ppp_SP <- as.SpatialPoints.ppp(zn1a)
summary(rev_ppp_SP)
rev_ppp_SPDF <- as.SpatialPointsDataFrame.ppp(zn1a)
summary(rev_ppp_SPDF)
rev_ppp_SGDF <- as.SpatialGridDataFrame.ppp(zn1a)
summary(rev_ppp_SGDF)
# as.owin(meuse.grid)
# --- will only work if as.owin is method...
mg_owin <- as(as(meuse.grid["ffreq"], "SpatialPixelsDataFrame"), "owin")
mg_owin
ho_sp <- SpatialPolygons(list(Polygons(list(Polygon(cbind(c(0,1,1,0,0),
  c(0,0,1,1,0))), Polygon(cbind(c(0.6,0.4,0.4,0.6,0.6),
  c(0.2,0.2,0.4,0.4,0.2)), hole=TRUE)), ID="ho")))
plot(ho_sp, col="red", pbg="pink")
ho <- as(ho_sp, "owin")
plot(ho)
pp <- runifpoint(500, win=ho)
plot(pp)
ho_orig <- owin(poly=list(list(x=c(0,1,1,0), y=c(0,0,1,1)),
  list(x=c(0.6,0.4,0.4,0.6), y=c(0.2,0.2,0.4,0.4))))
identical(ho, ho_orig)
mg_dist <- meuse.grid["dist"]
fullgrid(mg_dist) <- TRUE
image(mg_dist, axes=TRUE)
mg_im <- as(mg_dist, "im")
plot(mg_im)
mg2 <- as.SpatialGridDataFrame.im(mg_im)
image(mg2, axes=TRUE)
[Package maptools version 0.6-9 Index]


read/write to/from (ESRI) asciigrid format

Description

read/write to/from ESRI asciigrid format; a fuzz factor has been added to writeAsciiGrid to force cell resolution to equality if the difference is less than the square root of machine precision

Usage

readAsciiGrid(fname, as.image = FALSE, plot.image = FALSE,
 colname = basename(fname), proj4string = CRS(as.character(NA)),
 dec=options()$OutDec)
writeAsciiGrid(x, fname, attr = 1, na.value = -9999, dec=options()$OutDec, ...)

Arguments

fname | file name
as.image | logical; if TRUE, a list is returned, ready to be shown with the image command; if FALSE an object of class SpatialGridDataFrame-class is returned
plot.image | logical; if TRUE, an image of the map is plotted
colname | alternative name for data column if not file basename
proj4string | A CRS object setting the projection arguments of the Spatial Grid returned
dec | decimal point character. This should be a character string containing just one single-byte character — see note below.
x | object of class SpatialGridDataFrame
attr | attribute column; if missing, the first column is taken; a name or a column number may be given
na.value | numeric; value given to missing valued cells in the resulting map
... | arguments passed to write.table, which is used to write the numeric data

Value

readAsciiGrid returns the grid map read; either as an object of class SpatialGridDataFrame-class or, if as.image is TRUE, as list with components x, y and {z}.

Note

In ArcGIS 8, it was not in general necessary to set the dec argument; it is not necessary in a mixed environment with ArcView 3.2 (R writes and ArcView reads "."), but inter-operation with ArcGIS 9 requires care because the defaults used by ArcGIS seem to be misleading, and it may be necessary to override what appear to be platform defaults by setting the argument.

Author(s)

Edzer J. Pebesma, e.pebesma@geo.uu.nl

See Also

as.image.SpatialGridDataFrame, image

Examples

x <- readAsciiGrid(system.file("grids/test.ag", package="maptools")[1])
summary(x)
image(x)
xp <- as(x, "SpatialPixelsDataFrame")
abline(h=332000, lwd=3)
xpS <- xp[coordinates(xp)[,2] < 332000,]
summary(xpS)
xS <- as(xpS, "SpatialGridDataFrame")
summary(xS)
tmpfl <- paste(tempdir(), "testS.ag", sep="/")
writeAsciiGrid(xS, tmpfl)
axS <- readAsciiGrid(tmpfl)
opar <- par(mfrow=c(1,2))
image(xS, main="before export")
image(axS, main="after import")
par(opar)
unlink(tmpfl)


Check holes in Polygons objects

Description

The function checks holes in Polygons objects, using an intersection between a gpclib package gpc.poly object with one or more polygon contours and its bounding box to set the hole flag. The function will set single polygon contours to hole=FALSE, and if multiple polygon contours are holes, will set them TRUE.

Usage

checkPolygonsHoles(x)

Arguments

x An Polygons object as defined in package sp

Value

An Polygons object re-created from the input object.

Author(s)

Roger Bivand

Examples

library(sp)
library(gpclib)
nc1 <- readShapePoly(system.file("shapes/sids.shp", package="maptools")[1], proj4string=CRS("+proj=longlat +ellps=clrk66"))
pl <- getSpPpolygonsSlot(nc1)
sapply(getPolygonsPolygonsSlot(pl[[4]]), function(x) x@hole)
pl[[4]] <- Polygons(list(getPolygonsPolygonsSlot(pl[[4]])[[1]], Polygon(getPolygonCoordsSlot(getPolygonsPolygonsSlot(pl[[4]])[[2]]), hole=TRUE), getPolygonsPolygonsSlot(pl[[4]])[[3]]), getPolygonsIDSlot(pl[[4]]))
sapply(getPolygonsPolygonsSlot(pl[[4]]), function(x) x@hole)
pl_new <- lapply(pl, checkPolygonsHoles)
sapply(getPolygonsPolygonsSlot(pl_new[[4]]), function(x) x@hole)
srs <- getPolygonCoordsSlot(getPolygonsPolygonsSlot(pl[[1]])[[1]])
hle2 <- structure(c(-81.64093, -81.38380, -81.34165, -81.66833, -81.64093,
 36.57865, 36.57234, 36.47603, 36.47894, 36.57865), .Dim = as.integer(c(5, 2)))
hle3 <- structure(c(-81.47759, -81.39118, -81.38486, -81.46705, -81.47759,
 36.56289, 36.55659, 36.49907, 36.50380, 36.56289), .Dim = as.integer(c(5, 2)))
x <- Polygons(list(Polygon(srs), Polygon(hle2), Polygon(hle3)),
 ID=getPolygonsIDSlot(pl[[1]]))
sapply(getPolygonsPolygonsSlot(x), function(x) x@hole)
res <- checkPolygonsHoles(x)
sapply(getPolygonsPolygonsSlot(res), function(x) x@hole)
## Not run:
opar <- par(mfrow=c(1,2))
plot(SpatialPolygons(list(x)))
text(t(sapply(getPolygonsPolygonsSlot(x), getPolygonLabptSlot)),
 labels=sapply(getPolygonsPolygonsSlot(x), getPolygonHoleSlot), cex=0.6)
title(xlab="Hole slot values before checking")
plot(SpatialPolygons(list(res)))
text(t(sapply(getPolygonsPolygonsSlot(res), getPolygonLabptSlot)),
 labels=sapply(getPolygonsPolygonsSlot(res), getPolygonHoleSlot), cex=0.6)
title(xlab="Hole slot values after checking")
par(opar)
## End(Not run)


Converter functions to build SpatialLinesDataFrame objects

Description

These functions show how to build converters to SpatialLinesDataFrame objects: ArcObj2SLDF from the list returned by the get.arcdata function in the RArcInfo package; ContourLines2SLDF from the list returned by the contourLines function in the graphics package (here the data frame is just the contour levels, with one Lines object made up of at least one Line object per level); and MapGen2SL reads a file in "Mapgen" format into a SpatialLines object.

Usage

ArcObj2SLDF(arc, proj4string=CRS(as.character(NA)), IDs)
ContourLines2SLDF(cL, proj4string=CRS(as.character(NA)))
MapGen2SL(file, proj4string=CRS(as.character(NA)))

Arguments

arc | a list returned by the get.arcdata function in the RArcInfo package
IDs | vector of unique character identifiers; if not given, suitable defaults will be used, and the same values inserted as data slot row names
cL | a list returned by the contourLines function in the graphics package
proj4string | Object of class "CRS"; see CRS-class
file | filename of a file containing a Mapgen line data set

Value

A SpatialLinesDataFrame object

Note

Coastlines of varying resolution may be chosen online and downloaded in "Mapgen" text format from http://www.ngdc.noaa.gov/mgg/shorelines/shorelines.html, most conveniently using the interactive selection tool, but please note the 500,000 point limit on downloads, which is easy to exceed.

Author(s)

Roger Bivand; Edzer Pebesma

See Also

SpatialLines-class

Examples

#data(co37_d90_arc) # retrieved as:
# library(RArcInfo)
# fl <- "http://www.census.gov/geo/cob/bdy/co/co90e00/co37_d90_e00.zip"
# download.file(fl, "co37_d90_e00.zip")
# e00 <- zip.file.extract("co37_d90.e00", "co37_d90_e00.zip")
# e00toavc(e00, "ncar")
# arc <- get.arcdata(".", "ncar")
#res <- arcobj2SLDF(arc)
#plot(res)
#invisible(title(""))
res <- ContourLines2SLDF(contourLines(volcano))
plot(res, col=terrain.colors(nrow(as(res, "data.frame"))))
title("Volcano contours as SpatialLines")


Put dots in polygons

Description

Make point coordinates for a dot density map

Usage

dotsInPolys(pl, x, f = "random", offset, compatible = FALSE)

Arguments

pl | list of polygons of class polylist, or an object of class SpatialPolygons or SpatialPolygonsDataFrame
x | integer vector of counts of same length as pl for dots
f | type of sampling used to place points in polygons, either "random" or "regular"
offset | for regular sampling only: the offset (position) of the regular grid; if not set, c(0.5,0.5), that is the returned grid is not random
compatible | what to return, if TRUE a a list of matrices of point coordinates, one matrix for each member of pl, if false a SpatialPointsDataFrame with polygon ID values

Details

With f="random", the dots are placed in the polygon at random, f="regular" - in a grid pattern (number of dots not guaranteed to be the same as the count). When the polygon is made up of more than one part, the dots will be placed in proportion to the relative areas of the clockwise rings (anticlockwise are taken as holes). From maptools release 0.5-2, correction is made for holes in the placing of the dots, but depends on hole values being correctly set, which they often are not. The wrapper package spgpc may be used to check holes, see the dontrun section of the example.

Value

If compatible=TRUE, the function returns a list of matrices of point coordinates, one matrix for each member of pl. If x[i] is zero, the list element is NULL, and can be tested when plotting - see the examples. If compatible=FALSE (default), it returns a SpatialPointsDataFrame with polygon ID values as the only column in the data slot.

Note

Waller and Gotway (2004) Applied Spatial Statistics for Public Health Data (Wiley, Hoboken, NJ) explicitly warn that care is needed in plotting and interpreting dot density maps (pp. 81-83)

Author(s)

Roger Bivand Roger.Bivand@nhh.no

See Also

spsample

Examples

x <- read.shape(system.file("shapes/sids.shp", package="maptools")[1])
ncpolys <- Map2poly(x)
try1 <- dotsInPolys(ncpolys, as.integer(x$att.data$SID74), compatible=TRUE)
plot(ncpolys)
xx <- lapply(try1, function(x) {if (!is.null(x)) points(as.matrix(rbind(x)), pch=18, col="red")}) # added as.matrix(rbind(x)), thanks to Nicholas Grassly
try2 <- dotsInPolys(ncpolys, as.integer(x$att.data$SID74), f="regular", compatible=TRUE)
plot(ncpolys)
xx <- lapply(try2, function(x) {if (!is.null(x)) points(as.matrix(rbind(x)), pch=18, col="red")})
nc_SP <- readShapePoly(system.file("shapes/sids.shp", package="maptools")[1], proj4string=CRS("+proj=longlat  +ellps=clrk66"))
## Not run:
library(spgpc)
pls <- getSpPpolygonsSlot(nc_SP)
pls_new <- lapply(pls, checkPolygonsHoles)
nc_SP <- SpatialPolygonsDataFrame(SpatialPolygons(pls_new, proj4string=CRS(proj4string(nc_SP))), data=as(nc_SP, "data.frame"))
## End(Not run)
try1 <- dotsInPolys(nc_SP, as.integer(nc_SP$SID74))
plot(nc_SP, axes=TRUE)
plot(try1, add=TRUE, pch=18, col="red")
try2 <- dotsInPolys(nc_SP, as.integer(nc_SP$SID74), f="regular")
plot(nc_SP, axes=TRUE)
plot(try2, add=TRUE, pch=18, col="red")


Create SpatialGrid for PNG output to GE

Description

The function sets up metadata in the form of a SpatialGrid object for defining the size and placing of a PNG image overlay in Google Earth.

Usage

GE_SpatialGrid(obj, asp = NA, maxPixels = 600)

Arguments

obj | a Spatial* object
asp | if NA, will be set to the latitude corrected value
maxPixels | the maximum dimension of the output PNG

Details

The function is used together with kmlOverlay to wrap around the opening of a PNG graphics device, plotting code, and the closing of the device. The computed values take account of the adjustment of the actual data bounding box to an integer number of rows and columns in the image file.

The approach may be used as an alternative to writing PNG files from SpatialGrid and SpatialPixel objects in rgdal using writeGDAL, and to writing KML files using writeOGR for vector data objects. The output PNG files are likely to be very much smaller than large vector data KML files, and hinder the retrieval of exact positional information.

Value

returns an S3 object of class \dQuote{GE_SG} with components:

height | Integer raster height for png call
width | Integer raster width for png call
SG | a SpatialGrid object with the grid topology of the output PNG
asp | the aspect value used
xlim | xlim taken from SG
ylim | ylim taken from SG

Author(s)

Duncan Golicher, David Forrest and Roger Bivand

See Also

kmlOverlay

Examples

opt_exask <- options(example.ask=FALSE)
qk <- SpatialPointsDataFrame(quakes[, c(2:1)], quakes)
proj4string(qk) <- CRS("+proj=longlat")
tf <- tempfile()
SGqk <- GE_SpatialGrid(qk)
png(file=paste(tf, ".png", sep=""), width=SGqk$width, height=SGqk$height,
  bg="transparent")
par(mar=c(0,0,0,0), xaxs="i", yaxs="i")
plot(qk, xlim=SGqk$xlim, ylim=SGqk$ylim, setParUsrBB=TRUE)
dev.off()
kmlOverlay(SGqk, paste(tf, ".kml", sep=""), paste(tf, ".png", sep=""))
## Not run:
qk0 <- quakes
qk0$long <- ifelse(qk0$long <= 180, qk0$long, qk0$long-360)
qk0a <- SpatialPointsDataFrame(qk0[, c(2:1)], qk0)
proj4string(qk0a) <- CRS("+proj=longlat")
writeOGR(qk0a, paste(tf, "v.kml", sep=""), "Quakes", "KML")
system(paste("googleearth ", tf, ".kml", sep=""))
## End(Not run)
options(example.ask=opt_exask)


Polygon centroids

Description

return the centroids of a map of polygons

Usage

get.Pcent(theMap)

Arguments

theMap a Map object returned by read.shape()

Value

a matrix of centroids.

Author(s)

Nicholas J. Lewin-Koh, modified by Roger Bivand Roger.Bivand@nhh.no

Examples

x <- read.shape(system.file("shapes/sids.shp", package="maptools")[1])
get.Pcent(x)


Find azimuth for geographical coordinates

Description

The function finds azimuth values for geographical coordinates given as decimal degrees from the from coordinates to the to coordinate. In function trackAzimuth, the azimuth values are found between successive rows of the input coordinate matrix.

Usage

gzAzimuth(from, to, type = "snyder_sphere")
trackAzimuth(track, type="snyder_sphere")

Arguments

from | a two column matrix of geographical coordinates given as decimal degrees (longitude first)
track | a two column matrix of geographical coordinates given as decimal degrees (longitude first)
to | a one row, two column matrix or two element vector of geographical coordinates given as decimal degrees (longitude first)
type | default is "snyder_sphere", otherwise "abdali"; the results should be identical with slightly less trigonometry in "abdali"

Details

The azimuth is calculated on the sphere, using the formulae given by Snyder (1987, p. 30) and Abdali (1997, p. 17). The examples use data taken from Abdali (p. 17–18). There is a very interesting discussion of the centrality of azimuth-finding in the development of mathematics and mathematical geography in Abdali's paper. Among others, al-Khwarizmi was an important contributor. As Abdali puts it, "This is a veritable who's who of medieval science" (p. 3).

Value

values in decimal degrees - zero is North - of the azimuth from the from coordinates to the to coordinate.

Author(s)

Roger Bivand, with contributions by Sebastian Luque

References

Snyder JP (1987) Map projections - a working manual, USGS Professional Paper 1395; Abdali SK (1997) "The Correct Qibla", http://patriot.net/users/abdali/ftp/qibla.pdf

Examples

name <- c("Mecca", "Anchorage", "Washington")
long <- c(39.823333, -149.883333, -77.0166667)
lat <- c(21.423333, 61.2166667, 38.9)
x <- cbind(long, lat)
rownames(x) <- name
crib <- c(-9.098363, 56.575960)
r1 <- gzAzimuth(x[2:3,], x[1,])
r1
all.equal(r1, crib)
r2 <- gzAzimuth(x[2:3,], x[1,], type="abdali")
r2
all.equal(r2, crib)
trackAzimuth(x)


Hisaji Ono's lake/hole problem

Description

How to plot polygons with holes - holes are encoded by coordinates going anticlockwise, and overplotting is avoided by re-ordering the order in which polygons are plotted.

Usage

data(gpcholes)

Format

The format is class "polylist".

Details

"Date: Tue, 11 May 2004 12:54:20 +0900 From: Hisaji ONO To: r-help

I've tried to create a polygon with one hole by gpclib using following example script.

holepoly <- read.polyfile(system.file("poly-ex/hole-poly.txt", package="gpclib"), nohole = FALSE)
area.poly(holepoly)
plot(holepoly,poly.args=list(col="red",border="blue"))

And I noticed plot function couldn't draw polygons with holes correctly.

Does anyone know how to solve this situation?"

*(hole1pl has reversed the y component of polygon 1, to make its ring direction clockwise, hole2pl reverses the order of the two polygons in holepoly1@pts)*

Source

Data file included in "gpclib" package.

Examples

data(gpcholes)
plot(hole2pl, col="red", pbg="white", border="blue", forcefill=FALSE)
plot(hole1pl, col="red", pbg="white", border="blue", forcefill=FALSE)


Create and write KML file for PNG image overlay

Description

The function is used to create and write a KML file for a PNG image overlay for Google Earth.

Usage

kmlOverlay(obj, kmlfile = NULL, imagefile = NULL, name = "R image")

Arguments

obj | a \dQuote{GE_SG} object from GE_SpatialGrid
kmlfile | if not NULL the name of the kml file to be written
imagefile | the name of the PNG file containing the image - this should be either relative (same directory as kml file) or abosolute (fully qualified)
name | the name used to describe the image overlay in GE

Details

The function is used together with GE_SpatialGrid to wrap around the opening of a PNG graphics device, plotting code, and the closing of the device. The computed values take account of the adjustment of the actual data bounding box to an integer number of rows and columns in the image file.

The approach may be used as an alternative to writing PNG files from SpatialGrid and SpatialPixel objects in rgdal using writeGDAL, and to writing KML files using writeOGR for vector data objects. The output PNG files are likely to be very much smaller than large vector data KML files, and hinder the retrieval of exact positional information.

Value

x is a character vector containing the generated lines of the kml file

Author(s)

Duncan Golicher, David Forrest and Roger Bivand

See Also

GE_SpatialGrid

Examples

opt_exask <- options(example.ask=FALSE)
qk <- SpatialPointsDataFrame(quakes[, c(2:1)], quakes)
proj4string(qk) <- CRS("+proj=longlat")
tf <- tempfile()
SGqk <- GE_SpatialGrid(qk)
png(file=paste(tf, ".png", sep=""), width=SGqk$width, height=SGqk$height,
  bg="transparent")
par(mar=c(0,0,0,0), xaxs="i", yaxs="i")
plot(qk, xlim=SGqk$xlim, ylim=SGqk$ylim, setParUsrBB=TRUE)
dev.off()
kmlOverlay(SGqk, paste(tf, ".kml", sep=""), paste(tf, ".png", sep=""))
## Not run:
library(rgdal)
qk0 <- quakes
qk0$long <- ifelse(qk0$long <= 180, qk0$long, qk0$long-360)
qk0a <- SpatialPointsDataFrame(qk0[, c(2:1)], qk0)
proj4string(qk0a) <- CRS("+proj=longlat")
writeOGR(qk0a, paste(tf, "v.kml", sep=""), "Quakes", "KML")
system(paste("googleearth ", tf, ".kml", sep=""))
## End(Not run)
options(example.ask=opt_exask)


Create polygon lists and bounding boxes from imported shapefiles

Description

Map2poly() is a function to make imported GIS vector polygons into a "polylist" object from a "Map" object. Map2bbs() retrieves polygon bounding boxes; analogous functions for the "shapefiles" package are shape2poly() and shape2bbs() (thanks to Stéphane Dray for his contribution); convert.pl() serves to convert the deprecated "multipart" "polylist" format to an NA-separated format now used by plot() for polylist objects and poly2nb(). *2lines() and *2points() do the same as *2poly() for shapefiles of types 3=lines and 1=points; polygons are type=5.

Usage

Map2poly1(Map, region.id = NULL, raw=TRUE)
Map2poly(Map, region.id = NULL, quiet=TRUE)
Map2lines(Map)
Map2points(Map)
Map2bbs(Map)
MapShapeIds(Map)
shape2poly(shape, region.id = NULL)
shape2lines(shape)
shape2points(shape)
shape2bbs(shape)
convert.pl(pl)

Arguments

Map | object of class "Map", imported from a GIS polygon vector file
shape | list returned by read.shapefile(), imported from a GIS polygon vector file using the package "shapefiles"
region.id | character vector of region ids to be added to the neighbours list as attribute region.id
raw | default TRUE until next release, if TRUE, do not run sanity check for ring directions implying holes in no surrounding polygon
quiet | if TRUE, suppress reports of ring direction changes
pl | list of old-style multipart polygons of class polylist

Details

From release 0.4-12, Map2poly replaces the plotting order heuristic with plotting only in strict decreasing order of top-level (multi)polygon area. In multipolygons, the components are plotted in decreasing area order. The multpolygons themselves are assigned their area sum for plotting order. Ring directions are all set to clockwise - very many shapefiles have been observed with quite unrealistic ring orders, so respecting input ring orders seems worse than imposing uniformity. The previous plot order assignment mechanism is retained using Map2poly1. From release 0.4-1, Map2poly tries to determine a plotting order for multiple parts of a single polygon object, and for lists of polygon objects to avoid overplotting. It also tries to respect the ESRI-stipulated rule that lakes are anti-clockwise and islands are clockwise: http://shapelib.maptools.org/dl/shapefile.pdf: "The order of vertices or orientation for a ring indicates which side of the ring is the interior of the polygon. The neighborhood to the right of an observer walking along the ring in vertex order is the neighborhood inside the polygon. Vertices of rings defining holes in polygons are in a counterclockwise direction. Vertices for a single, ringed polygon are, therefore, always in clockwise order" (p.8). If raw=FALSE, an additional attempt is made to enforce the final stipulation, where the polygon's bounding box does not fall wholly inside the bounding box of another polygon. This may generate warnings indicating when this has taken place, indicating an unexpected ring direction in the imported shapefile.

Value

Map2poly, shape2poly, convert.pl return polylist objects, lists of polygon boundary coordinates (divided by NA where the polygon object is represented by more than one polygon) with many attributes; Map2lines, shape2lines return lineslist objects; Map2points, shape2points return (n x 2) matrices; and Map2bbs, shape2bbs return bounding box matrixes, c(x1, y1, x2, y2) with coordinates for the lower left corner and upper right corner.

Author(s)

Roger Bivand and Stéphane Dray, Roger.Bivand@nhh.no

See Also

read.shape, read.shapefile

Examples

try2 <- read.shape(system.file("shapes/columbus.shp", package="maptools")[1])
mappolys <- Map2poly(try2, as.character(try2$att.data$NEIGNO), quiet=FALSE)
plot(mappolys)
title(main="Polygons for Columbus OH from maptools package")
mappolys <- Map2poly1(try2, as.character(try2$att.data$NEIGNO), raw=FALSE)
plot(mappolys)
title(main="Polygons for Columbus OH from maptools package")
try3 <- read.shape(system.file("shapes/baltim.shp", package="maptools")[1])
baltpts <- Map2points(try3)
xylims <- attr(baltpts, "maplim")
plot(xylims$x, xylims$y, asp=1, type='n', xlab="", ylab="")
points(baltpts)
title(main="Baltimore points from maptools package")
try4 <- read.shape(system.file("shapes/fylk-val.shp", package="maptools")[1])
fylk.val <- Map2lines(try4)
xylims <- attr(fylk.val, "maplim")
plot(xylims$x, xylims$y, asp=1, type='n', xlab="", ylab="")
for (i in 1:length(fylk.val)) lines(fylk.val[[i]])
title(main="Norwegian river centrelines from maptools package")


Convert map objects to sp classes

Description

These functions may be used to convert map objects returned by the map function in the maps package to suitable objects defined in the sp package. In the examples below, arguments are shown for retrieving first polygons by name, then lines by window.

Usage

map2SpatialPolygons(map, IDs, proj4string = CRS(as.character(NA)))
map2SpatialLines(map, IDs=NULL, proj4string = CRS(as.character(NA)))
pruneMap(map, xlim=NULL, ylim=NULL)

Arguments

map | a map object defined in the maps package and returned by the map function
IDs | Unique character ID values for each output Polygons object; the input IDs can be an integer or character vector with duplicates, where the duplicates will be combined as a single output Polygons object
proj4string | Object of class "CRS"; holding a valid proj4 string
xlim,ylim | limits for pruning a map object - should only be used for lines, because polygons will not be closed

Value

map2SpatialPolygons returns a SpatialPolygons object and map2SpatialLines returns a SpatialLines object (objects defined in the sp package); pruneMap returns a modified map object defined in the maps package

Note

As the examples show, retrieval by name should be checked to see whether a window is not also needed: the "norway" polygons include "Norway:Bouvet Island", which is in the South Atlantic. Here, the IDs argument is set uniformly to "Norway" for all the component polygons, so that the output object contains a single Polygons object with multiple component Polygon objects. When retrieving by window, pruning may be needed on lines which are included because they begin within the window; interior=FALSE is used to remove country boundaries in this case.

Author(s)

Roger Bivand

See Also

map

Examples

library(maps)
nor_coast_poly <- map("world", "norway", fill=TRUE, col="transparent", plot=FALSE)
range(nor_coast_poly$x, na.rm=TRUE)
range(nor_coast_poly$y, na.rm=TRUE)
nor_coast_poly <- map("world", "norway", fill=TRUE, col="transparent", plot=FALSE, ylim=c(58,72))
nor_coast_poly$names
IDs <- sapply(strsplit(nor_coast_poly$names, ":"), function(x) x[1])
nor_coast_poly_sp <- map2SpatialPolygons(nor_coast_poly, IDs=IDs, proj4string=CRS("+proj=longlat +datum=wgs84"))
getSpPnParts(nor_coast_poly_sp)
plot(nor_coast_poly_sp, col="grey", axes=TRUE)
nor_coast_lines <- map("world", interior=FALSE, plot=FALSE, xlim=c(4,32), ylim=c(58,72))
plot(nor_coast_lines, type="l")
nor_coast_lines <- pruneMap(nor_coast_lines, xlim=c(4,32), ylim=c(58,72))
lines(nor_coast_lines, col="red")
nor_coast_lines_sp <- map2SpatialLines(nor_coast_lines, proj4string=CRS("+proj=longlat +datum=wgs84"))
plot(nor_coast_poly_sp, col="grey", axes=TRUE)
plot(nor_coast_lines_sp, col="blue", add=TRUE)


Report version information and changes

Description

Access version information and changes

Usage

maptools(changes = FALSE)

Arguments

changes If TRUE, reports changes, most recent last

Value

Invisibly returns the change list, with members version and changes

Author(s)

Roger Bivand

Examples

maptools(changes=TRUE)


Break polygons at meridian for recentering

Description

When recentering a world map, say to change an "Atlantic" view with longitude range -180 to 180, to a "Pacific" view, with longitude range 0 to 360, polygons crossed by the new offset, here 0/360, need to be clipped into left and right sub.polygons to avoid horizontal scratches across the map. The nowrapSpatialPolygons function performs this operation using polygon intersection, and nowrapRecenter recenters the output SpatialPolygons object.

Usage

nowrapRecenter(obj, offset = 0, eps = rep(.Machine$double.eps, 2))
nowrapSpatialPolygons(obj, offset = 0, eps=rep(.Machine$double.eps, 2))

Arguments

obj | A SpatialPolygons object
offset | offset from the Greenwich meridian
eps | vector of two (left and right) fuzz factors to retract the ring from the offset

Value

A SpatialPolygons object

Author(s)

Roger Bivand

See Also

recenter, nowrapSpatialLines

Examples

library(gpclib)
## Not run:
library(maps)
world <- map("world", fill=TRUE, col="transparent", plot=FALSE)
worldSpP <- map2SpatialPolygons(world, world$names, CRS("+proj=longlat"))
worldSpPr <- recenter(worldSpP)
plot(worldSpPr)
title("Pacific view without polygon splitting")
worldSpPnr <- nowrapRecenter(worldSpP)
plot(worldSpPnr)
title("Pacific view with polygon splitting")
## End(Not run)
crds <- matrix(c(-1, 1, 1, -1, 50, 50, 52, 52), ncol=2)
rcrds <- rbind(crds, crds[1,])
SR <- SpatialPolygons(list(Polygons(list(Polygon(rcrds)), ID="r1")), proj4string=CRS("+proj=longlat"))
bbox(SR)
SRr <- recenter(SR)
bbox(SRr)
SRnr <- nowrapRecenter(SR)
bbox(SRnr)


Making SpatialPolygons objects from RArcInfo input

Description

This function is used in making SpatialPolygons objects from RArcInfo input.

Usage

pal2SpatialPolygons(arc, pal, IDs, dropPoly1=TRUE, proj4string=CRS(as.character(NA)))

Arguments

IDs | Unique character ID values for each output Polygons object; the input IDs can be an integer or character vector with duplicates, where the duplicates will be combined as a single output Polygons object
proj4string | Object of class "CRS"; holding a valid proj4 string
arc | Object returned by get.arcdata
pal | Object returned by get.paldata
dropPoly1 | Should the first polygon in the AVC or e00 data be dropped; the first polygon is typically the compound boundary of the whole dataset, and can be detected by looking at the relative lengths of the list components in the second component of pal, which are the numbers of arcs making up the boundary of each polygon

Value

The functions return a SpatialPolygons object

Author(s)

Roger Bivand

Examples

nc1 <- readShapePoly(system.file("shapes/sids.shp", package="maptools")[1], ID="FIPS")
plot(nc1)
text(getSpPPolygonsLabptSlots(nc1), labels=getSpPPolygonsIDSlots(nc1), cex=0.6)
library(maps)
ncmap <- map("county", "north carolina", fill=TRUE, col="transparent", plot=FALSE)
IDs <- sapply(strsplit(ncmap$names, "[,:]"), function(x) x[2])
nc2 <- map2SpatialPolygons(ncmap, IDs)
plot(nc2)
text(getSpPPolygonsLabptSlots(nc2), labels=getSpPPolygonsIDSlots(nc2), cex=0.6)
library(RArcInfo)
td <- tempdir()
tmpcover <- paste(td, "nc", sep="/")
if (!file.exists(tmpcover)) e00toavc(system.file("share/co37_d90.e00", package="maptools")[1], tmpcover)
arc <- get.arcdata(td, "nc")
pal <- get.paldata(td, "nc")
pat <- get.tabledata(paste(td, "info", sep="/"), "NC.PAT")
sapply(pal[[2]], function(x) length(x[[1]]))
IDs <- paste(pat$ST[-1], pat$CO[-1], sep="")
nc3 <- pal2SpatialPolygons(arc, pal, IDs=IDs)
plot(nc3)
text(getSpPPolygonsLabptSlots(nc3), labels=getSpPPolygonsIDSlots(nc3), cex=0.6)


Plot a Map object (deprecated)

Description

This function is deprecated. It is difficult to maintain and there are several alternatives, either by converting Map objects to sp class objects or polylist etc. objects. (The function plots a map object directly from the imported shapefile data).

Usage

## S3 method for class 'Map':
plot(x, recs, auxvar = NULL, add = FALSE, fg = "gray", ol = "black",
  prbg = NULL, glyph = 16, color='red', type = "q", nclass = 5, ...)

Arguments

x | a Map object returned by read.shape()
recs | a vector with the shapes to plot, if missing defaults to all shapes; if fg is not equal to the length of recs, all fg are set to fg[1]
auxvar | if the Map has polygon shapes, a variable from which to derive polygon fill colours using the computed class intervals, must be same length as number of shapes
add | default FALSE, if TRUE add to existing plot
fg | foreground colour, can be used to pass through point and polygon colours, permitting the built-in class interval calculation to be avoided; if fg is not equal to the length of recs, all fg are set to fg[1]
ol | line/boundary colour; boundaries may be removed by setting to NA oe "transparent"
prbg | if not NULL, draw a rectangle around the plot in this colour
glyph | points plotted as this pch code
color | base colour for color.ramp(); default "red"
type | default "q" for quantile classes for color.ramp(), can also be "e" for equal intervals
nclass | number of classes for class intervals
... | passed through to plotting functions

Value

If the auxvar variable is not used, the function returns NULL, otherwise it returns the list constructed by maptools:::color.ramp() with components:

ramp | vector of colours
col.class | vector of classes
breaks | class intervals given as argument to cut

Author(s)

Nicholas J. Lewin-Koh, modified by Roger Bivand Roger.Bivand@nhh.no

See Also

read.shape, readShapePoly, readShapeLines, readShapePoints, getinfo.shape

Examples

## Not run:
x <- read.shape(system.file("shapes/sids.shp", package="maptools")[1])
plot(x)
nParts <- sapply(x$Shapes, function(x) attr(x, "nPart"))
table(nParts)
cols <- c("azure", "blue", "orange")
fgs <- cols[nParts]
plot(x, fg=fgs)
res <- plot(x, auxvar=x$att.data$BIR74)
str(res)
res <- plot(x, auxvar=x$att.data$BIR74, type="e")
str(res)
## End(Not run)


Plot polygons

Description

A helper function for plotting polygons in a global bounding box, attempts to handle overplotting internally, holes should be coded by coordinates ordered anti-clockwise. Earlier behaviour requires forcefill=TRUE; this is now default, but will cease to be so from next release. Changes caused by some imported shapefiles data leading to erroneously plotting lakes/islands, and R polygon() overplotting internal polygons. leglabs makes character strings from the same break points. The plot.polylist() function may be used as a generic S3 method. NOTE! plotpolys() DEPRECATED: use generic plot() for polylist objects function instead.

Usage

## S3 method for class 'polylist':
plot(x, col, border = par("fg"), add = FALSE, xlim=NULL,
 ylim=NULL, xlab = "", ylab = "", xpd = NULL, density = NULL, angle = 45,
 pbg=NULL, forcefill=TRUE, ...)
leglabs(vec, under="under", over="over", between="-")
plotpolys(pl, bb, col = NA, border = par("fg"), add = FALSE, xlim=NULL, ylim=NULL, ...)

Arguments

pl, x | list of polygons of class polylist
bb | matrix of polygon bounding boxes - columns are LL(x), LL(y), UR(x), UR(y); note that this argument may be omitted if the polygon list object has an "maplim" attribute, or if both xlim and ylim arguments are given
col | colurs to use for filling the polygons
border | the color to draw the border
add | add to existing plot
xlim, ylim | numeric of length 2, giving the x and y coordinates ranges
xlab, ylab | x and y axis labels, default no label
xpd | (where) should clipping take place?
density | the density of shading lines, in lines per inch
angle | the slope of shading lines, given as an angle in degrees (counter-clockwise)
pbg | colour to be used for hole and background fill, by default NULL (note that par("bg") may be "transparent")
forcefill | default TRUE until next release - fill anyway, if FALSE: believe ring direction as indication of holes/lakes
... | other arguments passed to plot to set the plot window - not passed to polygon
vec | vector of break values
under | character value for under
over | character value for over
between | character value for between

Author(s)

Roger Bivand Roger.Bivand@nhh.no

See Also

Map2poly, findInterval

Examples

try2 <- read.shape(system.file("shapes/columbus.shp", package="maptools")[1])
mappolys <- Map2poly(try2, as.character(try2$att.data$NEIGNO))
brks <- round(quantile(try2$att.data$CRIME, probs=seq(0,1,0.2)), digits=2)
colours <- c("salmon1", "salmon2", "red3", "brown", "black")
plot(mappolys, col=colours[findInterval(try2$att.data$CRIME, brks, all.inside=TRUE)], forcefill=FALSE)
legend(x=c(5.8, 7.1), y=c(13, 14.5), legend=leglabs(brks),
  fill=colours, bty="n")
invisible(title(main=paste("Columbus OH: residential burglaries and vehicle", "thefts per thousand households, 1980", sep="\n")))
try3 <- read.shape(system.file("shapes/sids.shp", package="maptools")[1])
mappolys <- Map2poly(try3, as.character(try3$att.data$FIPSNO))
cols <- c("pink", "orange", "red")
np <- sapply(mappolys, function(x) attr(x, "nParts"))
plot(mappolys, col=cols[np], forcefill=FALSE)
invisible(title(main="White: one part, orange: two part, red: three part polygons"))
plot(mappolys, density=c(10, 20, 30)[np], angle=c(-45, 0, 45)[np], axes=FALSE, forcefill=FALSE)


Virtual class "ppp"

ppp-class(maptools) R Documentation

Virtual class 'ppp'

Description

Virtual S4 class definition for S3 classes in the spatstat package to allow S4-style coercion to these classes

Objects from the Class

A virtual Class: No objects may be created from it.

Author(s)

Edzer J. Pebesma


Read shapefile into Map object

Description

Read shapefile into Map object; the file should be given including its ".shp" extension, and the function will reconstruct the names of the database (dbf) file and the index (shx) file from these.

Usage

read.shape(filen, dbf.data = TRUE, verbose=TRUE, repair=FALSE)
getinfo.shape(filen)
## S3 method for class 'shapehead':
print(x, ...)

Arguments

filen | name of file with *.shp extension
dbf.data | read DBF data together with shapes, default TRUE
verbose | default TRUE — report type of shapefile and number of shapes
repair | default FALSE: some shapefiles provided by Geolytics Inc. have values of object sizes stored in the *.shx index file that are eight bytes too large, leading the function to try to read past the end of file. If repair=TRUE, an attempt is made to repair the internal values, permitting such files to be read.
x | a shapehead list as returned by getinfo.shape
... | other arguments passed to print

Details

The function calls code from shapelib to read shapefiles, a file format used by ESRI GIS software among others

Value

read.shape() returns either a list of shapes of class ShapeList, or if dbf.data = TRUE a Map object with:

Shapes | a list of shapes of class ShapeList; both the individual shapes and the list have attributes
att.data | a data frame of data from the associated DBF file; note that the field names are adjusted to use in R using make.names(), and so will permit the underscore character from R release 1.9.0.

Author(s)

Nicholas J. Lewin-Koh, modified by Roger Bivand Roger.Bivand@nhh.no; shapelib by Frank Warmerdam

References

http://shapelib.maptools.org/

See Also

plot.Map, read.dbf

Examples

x <- read.shape(system.file("shapes/sids.shp", package="maptools")[1])
length(x$Shapes)
unlist(lapply(x$att.data, class))
str(getinfo.shape(system.file("shapes/fylk-val.shp", package="maptools")[1]))


GPSbabel read interface

Description

The function reads a data frame from an attached GPS using the external program gpsbabel. The columns of the data frame need to be identified by hand because different GPS order NMEA data in different ways, and the columns should be converted to the correct classes by hand. Once the specifics of a particular GPS are identified, and ways of cleaning erroneous locations are found, the conversion of the output data frame into a usable one may be automated.

Usage

readGPS(i = "garmin", f = "usb:")

Arguments

i | INTYPE: a supported file type, default "garmin"
f | INFILE: the appropriate device interface, default "usb:", on Windows for serial interfaces commonly "com4:" or similar

Details

The function just wraps: gpsbabel -i INTYPE -f INFILE -o tabsep -F - in system(), and reads the returned character vector of lines into a data frame. On some systems, INFILE may not be readable by ordinary users without extra configuration. The gpsbabel program must be present and on the user's PATH for the function to work. Typically, for a given GPS, the user will have to experiment first to find a set of data-cleaning tricks that work, but from then on they should be repeatable.

Value

A data frame of waypoint values

Author(s)

Patrick Giraudoux and Roger Bivand

References

http://www.gpsbabel.org

Examples

## Not run:
b1 <- readGPS(f="usb:")
str(b1)
b2 <- b1[1:172,]
wp0 <- b2[,c(2,3,4,8,9,19)]
str(wp0)
wp0$long <- wp0$V9
wp0$lat <- as.numeric(as.character(wp0$V8))
wp0$id <- as.character(wp0$V2)
wp0$alt <- as.numeric(substring(as.character(wp0$V19), 1, (nchar(as.character(wp0$V19))-1)))
wp0$time <- as.POSIXct(strptime(paste(as.character(wp0$V3), as.character(wp0$V4)), format="
str(wp0)
wp1 <- wp0[,-(1:6)]
str(wp1)
summary(wp1)
## End(Not run)


Read arc shape files into SpatialLinesDataFrame objects

Description

The readShapeLines function reads data from an arc/line shapefile into a SpatialLinesDataFrame object; the shapefile may be of type polygon, but for just plotting for example coastlines, a SpatialLines object is sufficient. The writeLinesShape function writes data from a SpatialLinesDataFrame object to a shapefile.

Usage

readShapeLines(fn, proj4string=CRS(as.character(NA)), verbose=FALSE, repair=FALSE)
writeLinesShape(x, fn, factor2char = TRUE)

Arguments

fn | shapefile layer name, without extension
proj4string | Object of class CRS; holding a valid proj4 string
verbose | default TRUE - report type of shapefile and number of shapes
repair | default FALSE: some shapefiles provided by Geolytics Inc. have values of object sizes stored in the *.shx index file that are eight bytes too large, leading the function to try to read past the end of file. If repair=TRUE, an attempt is made to repair the internal values, permitting such files to be read.
x | a SpatialLinesDataFrame object
factor2char | logical, default TRUE, convert factor columns to character

Details

The shpID values of the shapefile will be used as Lines ID values; when writing shapefiles, the object data slot rownames are added to the DBF file as column SL_ID.

Value

a SpatialLinesDataFrame object

Author(s)

Roger Bivand

Examples

xx <- readShapeLines(system.file("shapes/fylk-val.shp", package="maptools")[1], proj4string=CRS("+proj=utm +zone=33 +datum=WGS84"))
plot(xx, col="blue")
summary(xx)
xxx <- xx[xx$LENGTH > 30000,]
plot(xxx, col="red", add=TRUE)
tmpfl <- paste(tempdir(), "xxline", sep="/")
writeLinesShape(xxx, tmpfl)
getinfo.shape(paste(tmpfl, ".shp", sep=""))
axx <- readShapeLines(tmpfl, proj4string=CRS("+proj=utm +zone=33 +datum=WGS84"))
plot(xxx, col="black", lwd=4)
plot(axx, col="yellow", lwd=1, add=TRUE)
unlink(paste(tmpfl, ".*", sep=""))
xx <- readShapeLines(system.file("shapes/sids.shp", package="maptools")[1], proj4string=CRS("+proj=longlat +datum=NAD27"))
plot(xx, col="blue")


Read points shape files into SpatialPointsDataFrame objects

Description

The readShapePoints reads data from a points shapefile into a SpatialPointsDataFrame object. The writePointsShape function writes data from a SpatialPointsDataFrame object to a shapefile. Both reading and writing can be carried out for 2D and 3D point coordinates.

Usage

readShapePoints(fn, proj4string = CRS(as.character(NA)), verbose = FALSE, repair=FALSE)
writePointsShape(x, fn, factor2char = TRUE)

Arguments

fn | shapefile layer name, without extension
proj4string | Object of class CRS; holding a valid proj4 string
verbose | default TRUE - report type of shapefile and number of shapes
repair | default FALSE: some shapefiles provided by Geolytics Inc. have values of object sizes stored in the *.shx index file that are eight bytes too large, leading the function to try to read past the end of file. If repair=TRUE, an attempt is made to repair the internal values, permitting such files to be read.
x | a SpatialLinesDataFrame object
factor2char | logical, default TRUE, convert factor columns to character

Value

a SpatialPointsDataFrame object

Author(s)

Roger Bivand

Examples

library(maptools)
xx <- readShapePoints(system.file("shapes/baltim.shp", package="maptools")[1])
plot(xx)
summary(xx)
xxx <- xx[xx$PRICE < 40,]
tmpfl <- paste(tempdir(), "xxpts", sep="/")
writePointsShape(xxx, tmpfl)
getinfo.shape(paste(tmpfl, ".shp", sep=""))
axx <- readShapePoints(tmpfl)
plot(axx, col="red", add=TRUE)
unlink(paste(tmpfl, ".*", sep=""))
xx <- readShapePoints(system.file("shapes/pointZ.shp", package="maptools")[1])
dimensions(xx)
plot(xx)
summary(xx)


Read polygon shape files into SpatialPolygonsDataFrame objects

Description

The readShapePoly reads data from a polygon shapefile into a SpatialPolygonsDataFrame object. The writePolyShape function writes data from a SpatialPolygonsDataFrame object to a shapefile.

Usage

readShapePoly(fn, IDvar=NULL, proj4string=CRS(as.character(NA)), verbose=FALSE, repair=FALSE, force_ring=FALSE)
writePolyShape(x, fn, factor2char = TRUE)

Arguments

fn | shapefile layer name, without extension
IDvar | a character string: the name of a column in the shapefile DBF containing the ID values of the shapes - the values will be converted to a character vector
proj4string | Object of class CRS; holding a valid proj4 string
verbose | default TRUE - report type of shapefile and number of shapes
repair | default FALSE: some shapefiles provided by Geolytics Inc. have values of object sizes stored in the *.shx index file that are eight bytes too large, leading the function to try to read past the end of file. If repair=TRUE, an attempt is made to repair the internal values, permitting such files to be read.
force_ring | if TRUE, close unclosed input rings
x | a SpatialPolygonsDataFrame object
factor2char | logical, default TRUE, convert factor columns to character

Details

If no IDvar argument is given, the shpID values of the shapefile will be used as Polygons ID values; when writing shapefiles, the object data slot rownames are added to the DBF file as column SP_ID.

Value

a SpatialPolygonsDataFrame object

Author(s)

Roger Bivand

Examples

library(maptools)
xx <- readShapePoly(system.file("shapes/sids.shp", package="maptools")[1], IDvar="FIPSNO", proj4string=CRS("+proj=longlat +ellps=clrk66"))
plot(xx, border="blue", axes=TRUE, las=1)
text(getSpPPolygonsLabptSlots(xx), labels=getSpPPolygonsIDSlots(xx), cex=0.6)
as(xx, "data.frame")[1:5, 1:6]
xxx <- xx[xx$SID74 < 2,]
plot(xxx, border="red", add=TRUE)
tmpfl <- paste(tempdir(), "xxline", sep="/")
writePolyShape(xxx, tmpfl)
getinfo.shape(paste(tmpfl, ".shp", sep=""))
axx <- readShapePoly(tmpfl, proj4string=CRS("+proj=longlat +ellps=clrk66"))
plot(xxx, border="black", lwd=4)
plot(axx, border="yellow", lwd=1, add=TRUE)
unlink(paste(tmpfl, ".*", sep=""))


Read GSHHS polygons into SpatialPolygons object

Description

Read GSHHS polygons into SpatialPolygons object for a chosen region, using binary shorelines from Global Self-consistant Hierarchical High-resolution Shorelines. The data are provided in integer form as millionths of decimal degrees,

Usage

Rgshhs(fn, xlim = NULL, ylim = NULL, level = 4, minarea = 0, shift=FALSE,
verbose = TRUE)

Arguments

fn | filename or full path to GSHHS file to be read
xlim | longitude limits within 0-360 in most cases, negative longitudes are also found east of the Atlantic, but the Americas are recorded as positive values
ylim | latitude limits
level | maximum GSHHS level to include, defaults to 4 (everything), setting 1 will only retrieve land, no lakes
minarea | minimum area in square km to retrieve, default 0
shift | default FALSE, can be used to shift longitudes > 180 degrees to below zero, beware of artefacts involving unhandled polygon splitting at 180 degrees
verbose | default TRUE, print progress reports

Details

The package is distributed with the coarse version of the shoreline data, and much more detailed versions may be downloaded from the referenced websites. The data is of high quality, matching the accuracy of SRTM shorelines for the full dataset (but not for inland waterbodies). In general, users will construct study region SpatialPolygons objects, which can then be exported (for example as a shapefile), or used in other R packages (such as PBSmapping). The largest land polygons take considerable time to clip to the study region, certainly many minutes for an extract from the full resolution data file including Eurasia (with Africa) or the Americas. For this reason, do not give up if nothing seems to be happening after the (verbose) message: "Rgshhs: clipping <m> of <n> polygons ..." appears. Clipping the largest polygons in full resolution also needs a good deal of memory

Value

a list with the following components:

polydata | data from the headers of the selected GSHHS polygons
belongs | a matrix showing which polygon belongs to (is included in) which polygon, going from the highest level among the selected polygons down to 1 (land); levels are: 1 land, 2 lake, 3 island_in_lake, 4 pond_in_island_in_lake.
new_belongs | a ragged list of polygon inclusion used for making SP
SP | a SpatialPolygons object; this is the principal output object, and will become the only output object as the package matures

Note

A number of steps are taken in this implementation that are unexpected, print messages, and so require explanation. Following the extraction of polygons intersecting the required region, a check is made to see if Antarctica is present. If it is, a new southern border is imposed at the southern ylim value or -90 if no ylim value is given. When clipping polygons seeming to intersect the required region boundary, it can happen that no polygon is left within the region (for example when the boundaries are overlaid, but also because the min/max polygon values in the header may not agree with the polygon itself (one case observed for a lake west of Groningen). The function then reports a null polygon. Another problem occurs when closed polygons are cut up during the finding of intersections between polygons and the required region boundary. The code in gpclib does not close them, so they are closed later, and the closure noted.

Please also note that use of gpclib is not limited in any way; the licence limitations only apply to the compilation of GPC C code into commercial software. This can be verified by checking that exactly the same GPC code is included in the GPL'ed PBSmapping package on CRAN.

Author(s)

Roger Bivand

References

http://www.soest.hawaii.edu/wessel/gshhs/gshhs.html, http://www.ngdc.noaa.gov/mgg/shorelines/gshhs.html; Wessel, P., and W. H. F. Smith, A Global Self-consistent, Hierarchical, High-resolution Shoreline Database, J. Geophys. Res., 101, 8741-8743, 1996.

Examples

gshhs.c.b <- system.file("share/gshhs_c.b", package="maptools")
NZx <- c(160,180)
NZy <- c(-50,-30)
NZ <- Rgshhs(gshhs.c.b, xlim=NZx, ylim=NZy)
plot(NZ$SP, col="khaki", pbg="azure2", xlim=NZx, ylim=NZy, xaxs="i", yaxs="i", axes=TRUE)
GLx <- c(265,285)
GLy <- c(40,50)
GL <- Rgshhs(gshhs.c.b, xlim=GLx, ylim=GLy)
plot(GL$SP, col="khaki", pbg="azure2", xlim=GLx, ylim=GLy, xaxs="i", yaxs="i", axes=TRUE)


write map data for Mondrian

Description

The function outputs a SpatialPolygonsDataFrame object to be used by Mondrian

Usage

sp2Mondrian(SP, file)

Arguments

SP | a SpatialPolygonsDataFrame object
file | file where output is written

Note

At this release, the function writes out a text file with both data and polygon(s) identified as belonging to each row of data.

Author(s)

Patrick Hausmann and Roger Bivand

References

http://www.theusrus.de/Mondrian/index.html

Examples

## Not run:
xx <- readShapePoly(system.file("shapes/columbus.shp", package="maptools")[1])
sp2Mondrian(xx, file="columbus1.txt")
xx <- readShapePoly(system.file("shapes/sids.shp", package="maptools")[1]))
sp2Mondrian(xx, file="sids1.txt")
## End(Not run)


Convert SpatialPolygons object for Stata tmap command

Description

The function converts a SpatialPolygons object for use with the Stata tmap command, by creating a data frame with the required columns.

Usage

sp2tmap(SP)

Arguments

SP a SpatialPolygons object

Value

a data frame with three columns:

_ID | an integer vector of polygon identifiers in numeric order
_X | numeric x coordinate
_Y | numeric y coordinate

and an ID_n attribute with the named polygon identifiers

Author(s)

Roger Bivand

References

http://www.stata.com/search.cgi?query=tmap

See Also

write.dta

Examples

## Not run:
xx <- readShapePoly(system.file("shapes/sids.shp", package="maptools")[1], IDvar="FIPSNO", proj4string=CRS("+proj=longlat +ellps=clrk66"))
plot(xx, border="blue", axes=TRUE, las=1)
tmapdf <- sp2tmap(as(xx, "SpatialPolygons"))
write.dta(tmapdf, file="NCmap.dta", version=7)
NCdf <- as(xx, "data.frame")
NCdf$ID_n <- attr(tmapdf, "ID_names")
write.dta(NCdf, file="NC.dta", version=7)
## End(Not run)


Export SpatialPolygons object as S-Plus map for WinBUGS

Description

The function exports an sp SpatialPolygons object into a S-Plus map format to be import by WinBUGS.

Usage

sp2WB(map, filename, Xscale = 1, Yscale = Xscale, plotorder = FALSE)

Arguments

map | a SpatialPolygons object
filename | file where output is written
Xscale, Yscale | scales to be written in the output file
plotorder | default=FALSE, if TRUE, export polygons in plotting order

Author(s)

Virgilio Gómez Rubio, partly derived from earlier code by Thomas Jagger

References

http://www.mrc-bsu.cam.ac.uk/bugs/winbugs/geobugs12manual.pdf

Examples

## Not run:
xx <- readShapePoly(system.file("shapes/sids.shp", package="maptools")[1], IDvar="FIPSNO", proj4string=CRS("+proj=longlat +ellps=clrk66"))
plot(xx, border="blue", axes=TRUE, las=1)
sp2WB(as(xx, "SpatialPolygons"), filename="NC.map")
x <- readAsciiGrid(system.file("grids/test.ag", package="maptools")[1])
xp <- as(x, "SpatialPixelsDataFrame")
pp <- as.SpatialPolygons.SpatialPixels(xp)
sp2WB(pp, filename="test.map")
## End(Not run)


Convert sp line and polygon objects to PBSmapping PolySet objects

Description

Functions SpatialLines2PolySet and SpatialPolygons2PolySet convert objects of sp classes to PolySet class objects as defined in the PBSmapping package.

Usage

SpatialLines2PolySet(SL)
SpatialPolygons2PolySet(SpP)

Arguments

SL | a SpatialLines object as defined in the sp package
SpP | a SpatialPolygons object as defined in the sp package

Value

PolySet objects as defined in the PBSmapping package

Author(s)

Roger Bivand

See Also

as.PolySet, MapGen2SL

Examples

library(PBSmapping)
library(maps)
nor_coast_lines <- map("world", interior=FALSE, plot=FALSE, xlim=c(4,32), ylim=c(58,72))
nor_coast_lines <- pruneMap(nor_coast_lines, xlim=c(4,32), ylim=c(58,72))
nor_coast_lines_sp <- map2SpatialLines(nor_coast_lines, proj4string=CRS("+proj=longlat +datum=wgs84"))
nor_coast_lines_PS <- SpatialLines2PolySet(nor_coast_lines_sp)
summary(nor_coast_lines_PS)
plotLines(nor_coast_lines_PS)
nor_coast_poly <- map("world", "norway", fill=TRUE, col="transparent", plot=FALSE, ylim=c(58,72))
IDs <- sapply(strsplit(nor_coast_poly$names, ":"), function(x) x[1])
nor_coast_poly_sp <- map2SpatialPolygons(nor_coast_poly, IDs=IDs, proj4string=CRS("+proj=longlat +datum=wgs84"))
nor_coast_poly_PS <- SpatialPolygons2PolySet(nor_coast_poly_sp)
summary(nor_coast_poly_PS)
plotPolys(nor_coast_poly_PS)


cbind for spatial objects

Description

spCbind provides cbind-like methods for Spatial*DataFrame objects in addition to the $, [<- and [[<- methods already available.

Methods

obj = "SpatialPointsDataFrame", x = "data.frame"
cbind a data frame to the data slot of a SpatialPointsDataFrame object
obj = "SpatialPointsDataFrame", x = "vector"
cbind a vector to the data slot of a SpatialPointsDataFrame object
obj = "SpatialLinesDataFrame", x = "data.frame"
cbind a data frame to the data slot of a SpatialLinesDataFrame object; the data frame argument must have row names set to the Lines ID values, and should be re-ordered first by matching against a shared key column
obj = "SpatialLinesDataFrame", x = "vector"
cbind a vector to the data slot of a SpatialLinesDataFrame object
obj = "SpatialPolygonsDataFrame", x = "data.frame"
cbind a data frame to the data slot of a SpatialPolygonsDataFrame object; the data frame argument must have row names set to the Polygons ID values, and should be re-ordered first by matching against a shared key column
obj = "SpatialPolygonsDataFrame", x = "vector"
cbind a vector to the data slot of a SpatialPolygonsDataFrame object

Author(s)

Roger Bivand

See Also

spChFIDs-methods, spRbind-methods

Examples

xx <- readShapePoly(system.file("shapes/sids.shp", package="maptools")[1],
  IDvar="FIPSNO", proj4string=CRS("+proj=longlat +ellps=clrk66"))
xtra <- read.dbf(system.file("share/nc_xtra.dbf", package="maptools")[1])
o <- match(xx$CNTY_ID, xtra$CNTY_ID)
xtra1 <- xtra[o,]
row.names(xtra1) <- xx$FIPSNO
xx1 <- spCbind(xx, xtra1)
names(xx1)
identical(xx1$CNTY_ID, xx1$CNTY_ID.1)


change feature IDs in spatial objects

Description

When the feature IDs need to be changed in SpatialLines* or SpatialPolygons* objects, these methods may be used. The new IDs should be a character vector of unique IDs of the correct length.

Methods

obj = "SpatialLines", x = "character"
replace IDs in a SpatialLines object
obj = "SpatialLinesDataFrame", x = "character"
replace IDs in a SpatialLinesDataFrame object
obj = "SpatialPolygons", x = "character"
replace IDs in a SpatialPolygons object
obj = "SpatialPolygonsDataFrame", x = "character"
replace IDs in a SpatialPolygonsDataFrame object

Note

It is usually sensible to keep a copy of the original feature IDs in the object, but this should be done by the user.

Author(s)

Roger Bivand

See Also

spCbind-methods, spRbind-methods

Examples

xx <- readShapePoly(system.file("shapes/sids.shp", package="maptools")[1],
  IDvar="FIPSNO", proj4string=CRS("+proj=longlat +ellps=clrk66"))
row.names(as(xx, "data.frame"))
xx1 <- spChFIDs(xx, as.character(xx$CNTY_ID))
row.names(as(xx1, "data.frame"))


rbind for spatial objects

Description

spRbind provides rbind-like methods for Spatial*DataFrame objects

Methods

obj = "SpatialPoints", x = "SpatialPoints"
rbind two SpatialPoints objects
obj = "SpatialPointsDataFrame", x = "SpatialPointsDataFrame"
rbind two SpatialPointsDataFrame objects
obj = "SpatialLines", x = "SpatialLines"
rbind two SpatialLines objects
obj = "SpatialLinesDataFrame", x = "SpatialLinesDataFrame"
rbind two SpatialLinesDataFrame objects
obj = "SpatialPolygons", x = "SpatialPolygons"
rbind two SpatialPolygons objects
obj = "SpatialPolygonsDataFrame", x = "SpatialPolygonsDataFrame"
rbind two SpatialPolygonsDataFrame objects

Author(s)

Roger Bivand

See Also

spChFIDs-methods, spCbind-methods

Examples

xx <- readShapePoly(system.file("shapes/sids.shp", package="maptools")[1],
  IDvar="FIPSNO", proj4string=CRS("+proj=longlat +ellps=clrk66"))
summary(xx)
xx$FIPSNO
xx1 <- xx[xx$CNTY_ID < 1982,]
xx2 <- xx[xx$CNTY_ID >= 1982,]
xx3 <- spRbind(xx2, xx1)
summary(xx3)
xx3$FIPSNO


Subset polygon list objects

Description

The function subsets a polygon list object, also subsetting region ID attributes and also map limits if required.

Usage

## S3 method for class 'polylist':
subset(x, subset, fit.bbox = TRUE, ...)

Arguments

x | a polylist object
subset | a logical vector valued TRUE if the element is to be retained
fit.bbox | if TRUE (default) modifies the maplim attribute to bound the subset
... | other arguments passed through

Value

returns a polylist object, lists of polygon boundary coordinates (divided by NA where the polygon object is represented by more than one polygon);

Author(s)

Roger Bivand Roger.Bivand@nhh.no

Examples

nc <- read.shape(system.file("shapes/sids.shp", package = "maptools")[1])
mappolys <- Map2poly(nc, as.character(nc$att.data$FIPSNO))
submap <- subset(mappolys, nc$att.data$SID74 > 0)
plot(mappolys, col="orange")
plot(submap, add=TRUE, col="lightpink", forcefill=TRUE)


Methods for sun ephemerides calculations

Description

Functions for calculating sunrise, sunset, and times of dawn and dusk, with flexibility for the various formal definitions. They use algorithms provided by the National Oceanic & Atmospheric Administration (NOAA).

Usage

## S4 method for signature 'SpatialPoints, POSIXct':
crepuscule(crds, dateTime, solarDep, direction=c("dawn", "dusk"),
           POSIXct.out=FALSE)
## S4 method for signature 'matrix, POSIXct':
crepuscule(crds, dateTime,
           proj4string=CRS("+proj=longlat +datum=WGS84"), solarDep,
           direction=c("dawn", "dusk"), POSIXct.out=FALSE)
## S4 method for signature 'SpatialPoints, POSIXct':
sunriset(crds, dateTime, direction=c("sunrise", "sunset"),
         POSIXct.out=FALSE)
## S4 method for signature 'matrix, POSIXct':
sunriset(crds, dateTime,
         proj4string=CRS("+proj=longlat +datum=WGS84"),
         direction=c("sunrise", "sunset"), POSIXct.out=FALSE)
## S4 method for signature 'SpatialPoints, POSIXct':
solarnoon(crds, dateTime, POSIXct.out=FALSE)
## S4 method for signature 'matrix, POSIXct':
solarnoon(crds, dateTime, proj4string=CRS("+proj=longlat +datum=WGS84"),
          POSIXct.out=FALSE)
## S4 method for signature 'SpatialPoints, POSIXct':
solarpos(crds, dateTime, ...)
## S4 method for signature 'matrix, POSIXct':
solarpos(crds, dateTime, proj4string=CRS("+proj=longlat +datum=WGS84"), ...)

Arguments

crds | a SpatialPoints or matrix object, containing x and y coordinates (in that order).
dateTime | a POSIXct object with the date and time associated to calculate ephemerides for points given in crds.
solarDep | numeric vector with the angle of the sun below the horizon in degrees.
direction | one of "dawn", "dusk", "sunrise", or "sunset", indicating which ephemerides should be calculated.
POSIXct.out | logical indicating whether POSIXct output should be included.
proj4string | string with valid projection string describing the projection of data in crds.
... | other arguments passed through.

Details

NOAA used the reference below to develop their Sunrise/Sunset http://www.srrb.noaa.gov/highlights/sunrise/sunrise.html and Solar Position http://www.srrb.noaa.gov/highlights/sunrise/azel.html Calculators. The algorithms include corrections for atmospheric refraction effects.

Input can consist of one location and at least one POSIXct times, or one POSIXct time and at least one location. solarDep is recycled as needed.

Do not use the daylight savings time zone string for supplying dateTime, as many OS will not be able to properly set it to standard time when needed.

Value

crepuscule, sunriset, and solarnoon return a numeric vector with the time of day at which the event occurs, expressed as a fraction, if POSIXct.out is FALSE; otherwise they return a data frame with both the fraction and the corresponding POSIXct date and time. solarpos returns a matrix with the solar azimuth (in degrees from North), and elevation.

Warning

Compared to NOAA's original Javascript code, the sunrise and sunset estimates from this translation may differ by +/- 1 minute, based on tests using selected locations spanning the globe. This translation does not include calculation of prior or next sunrises/sunsets for locations above the Arctic Circle or below the Antarctic Circle.

Note

NOAA notes that “for latitudes greater than 72 degrees N and S, calculations are accurate to within 10 minutes. For latitudes less than +/- 72 degrees accuracy is approximately one minute.”

Author(s)

Sebastian P. Luque spluque@gmail.com, translated from Greg Pelletier's gpel461@ecy.wa.gov VBA code (available from http://www.ecy.wa.gov/programs/eap/models.html), who in turn translated it from original Javascript code by NOAA (see Details). Roger Bivand roger.bivand@nhh.no adapted the code to work with sp classes.

References

Meeus, J. (1991) Astronomical Algorithms. Willmann-Bell, Inc.

Examples

## Location of Helsinki, Finland, in decimal degrees,
## as listed in NOAA's website
hels <- matrix(c(24.97, 60.17), nrow=1)
Hels <- SpatialPoints(hels, proj4string=CRS("+proj=longlat +datum=WGS84"))
d041224 <- as.POSIXct("2004-12-24", tz="EET")
## Astronomical dawn
crepuscule(hels, d041224, solarDep=18, direction="dawn", POSIXct.out=TRUE)
crepuscule(Hels, d041224, solarDep=18, direction="dawn", POSIXct.out=TRUE)
## Nautical dawn
crepuscule(hels, d041224, solarDep=12, direction="dawn", POSIXct.out=TRUE)
crepuscule(Hels, d041224, solarDep=12, direction="dawn", POSIXct.out=TRUE)
## Civil dawn
crepuscule(hels, d041224, solarDep=6, direction="dawn", POSIXct.out=TRUE)
crepuscule(Hels, d041224, solarDep=6, direction="dawn", POSIXct.out=TRUE)
solarnoon(hels, d041224, POSIXct.out=TRUE)
solarnoon(Hels, d041224, POSIXct.out=TRUE)
solarpos(hels, as.POSIXct(Sys.time(), tz="EET"))
solarpos(Hels, as.POSIXct(Sys.time(), tz="EET"))
sunriset(hels, d041224, direction="sunrise", POSIXct.out=TRUE)
sunriset(Hels, d041224, direction="sunrise", POSIXct.out=TRUE)
## Using a sequence of dates
Hels_seq <- seq(from=d041224, length.out=365, by="days")
up <- sunriset(Hels, Hels_seq, direction="sunrise", POSIXct.out=TRUE)
down <- sunriset(Hels, Hels_seq, direction="sunset", POSIXct.out=TRUE)
day_length <- down$time - up$time
plot(Hels_seq, day_length, type="l")

## Using a grid of spatial points for the same point in time
grd <- GridTopology(c(-179,-89), c(1,1), c(359,179))
SP <- SpatialPoints(coordinates(grd),
                    proj4string=CRS("+proj=longlat +datum=WGS84"))
wint <- as.POSIXct("2004-12-21", tz="GMT")
win <- crepuscule(SP, wint, solarDep=6, direction="dawn")
SPDF <- SpatialGridDataFrame(grd,
 proj4string=CRS("+proj=longlat +datum=WGS84"),
 data=data.frame(winter=win))
image(SPDF, axes=TRUE, col=cm.colors(40))


Place grids of points over polygons

Description

Place grids of points over polygons with chosen density and/or symbols (suggested by Michael Wolf).

Usage

symbolsInPolys(pl, dens, symb = "+", compatible = FALSE)

Arguments

pl | list of polygons of class polylist
dens | number of symbol plotting points per unit area; either a single numerical value for all polygons, or a numeric vector the same length as pl with values for each polygon
symb | plotting symbol; either a single value for all polygons, or a vector the same length as pl with values for each polygon
compatible | what to return, if TRUE a a list of matrices of point coordinates, one matrix for each member of pl, with a symb attribute, if false a SpatialPointsDataFrame with a symb column

Details

The dots are placed in a grid pattern with the number of points per polygon being polygon area times density (number of dots not guaranteed to be the same as the count). When the polygon is made up of more than one part, the dots will be placed in proportion to the relative areas of the clockwise rings (anticlockwise are taken as holes). From maptools release 0.5-2, correction is made for holes in the placing of the dots, but depends on hole values being correctly set, which they often are not. The wrapper package spgpc may be used to check holes, see the dontrun section of the example.

Value

The function returns a list of matrices of point coordinates, one matrix for each member of pl; each matrix has a symb attribute that can be used for setting the pch argument for plotting. If the count of points for the given density and polygon area is zero, the list element is NULL, and can be tested when plotting - see the examples.

Note

Extension to plot pixmaps at the plotting points using addlogo() from the pixmap package is left as an exercise for the user.

Author(s)

Roger Bivand Roger.Bivand@nhh.no

See Also

spsample

Examples

x <- read.shape(system.file("shapes/sids.shp", package="maptools")[1])
ncpolys <- Map2poly(x)
np <- sapply(ncpolys, function(x) attr(x, "nPart"))
symbs <- c("-", "+", "x")
try1 <- symbolsInPolys(ncpolys, 100, symb=symbs[np], compatible=TRUE)
plot(ncpolys)
xx <- lapply(try1, function(x) {if (!is.null(x)) points(x, pch=attr(x, "symb"))})
nc_SP <- readShapePoly(system.file("shapes/sids.shp", package="maptools")[1], proj4string=CRS("+proj=longlat +ellps=clrk66"))
## Not run:
library(spgpc)
pls <- getSpPpolygonsSlot(nc_SP)
pls_new <- lapply(pls, checkPolygonsHoles)
nc_SP <- SpatialPolygonsDataFrame(SpatialPolygons(pls_new, proj4string=CRS(proj4string(nc_SP))), data=as(nc_SP, "data.frame"))
## End(Not run)
symbs <- c("-", "+", "x")
np <- getSpPnParts(nc_SP)
try1 <- symbolsInPolys(nc_SP, 100, symb=symbs[np])
plot(nc_SP, axes=TRUE)
plot(try1, add=TRUE, pch=as.character(try1$symb))


Aggregate Polygons in a SpatialPolygons object

Description

The function aggregates Polygons in a SpatialPolygons object, according to the IDs vector specifying which input Polygons belong to which output Polygons; internal boundaries are dissolved using the gpclib package union() function.

Usage

unionSpatialPolygons(SpP, IDs, threshold=NULL)

Arguments

SpP | A SpatialPolygons object as defined in package sp
IDs | A vector defining the output Polygons objects, equal in length to the length of the polygons slot of SpRs; it may be character, integer, or factor (try table(factor(IDs)) for a sanity check)
threshold | if not NULL, an area measure below which slivers will be discarded (some polygons have non-identical boundaries, for instance along rivers, generating slivers on union which are artefacts, not real sub-polygons)

Value

Returns an aggregated SpatialPolygons object named with the aggregated IDs values in their sorting order; see getSpPPolygonsIDSlots() of the output object to view the order.

Author(s)

Roger Bivand

Examples

library(sp)
library(gpclib)
nc1 <- readShapePoly(system.file("shapes/sids.shp", package="maptools")[1], proj4string=CRS("+proj=longlat +datum=NAD27"))
lps <- getSpPPolygonsLabptSlots(nc1)
ID <- cut(lps[,1], quantile(lps[,1]), include.lowest=TRUE)
reg4 <- unionSpatialPolygons(nc1, ID)
getSpPPolygonsIDSlots(reg4)


Write a arc-type shapefile

Description

The function calls code from shapelib to write an arc-type shapefile from a list of matrices of line coordinates with no NAs.

Usage

write.linelistShape(linelist, df, file, factor2char = TRUE, strictFilename=FALSE)

Arguments

linelist | a list of matrices of line coordinates
df | a data frame object
file | a file name of maximum 8 chararacters, numbers or the underscore symbol to be written, omitting the extensions *.shp, *.shx and *.dbf, which are added in the function
factor2char | logical, default TRUE, convert factor columns to character
strictFilename | if TRUE, impose file basename length limit of 8 characters

Details

The function calls code from shapelib to write an arc-type shapefile (both the geometry file with a *.shp extension, the index file with a *.shx extension, and the database file with a *.dbf extension - see write.dbf for details of the data frame export within this function.

Value

no return value.

Note

From maptools 0.4-10, this function is placed in the user-visible namespace on a trial basis, and reports of any malfunction should be sent to the package maintainer, Roger Bivand Roger.Bivand@nhh.no. It is likely that this function and its arguments will be changed.

Author(s)

Nicholas J. Lewin-Koh, modified by Roger Bivand; shapelib by Frank Warmerdam

References

http://shapelib.maptools.org/

See Also

write.pointShape, write.dbf

Examples

x <- 10 * 1:nrow(volcano)
y <- 10 * 1:ncol(volcano)
line.list <- contourLines(x, y, volcano)
vol.levels <- data.frame(alt=sapply(line.list, function(x) x[[1]]))
vol.ll <- lapply(line.list, function(x) cbind(x$x, x$y))
tmpshp <- paste(tempdir(), "volcano", sep="/")
write.linelistShape(vol.ll, vol.levels, file=tmpshp)
try1 <- readShapeLines(tmpshp)
plot(try1)


Write a point-type shapefile

Description

The function calls code from shapelib to write a point-type shapefile.

Usage

write.pointShape(coordinates, df, file, factor2char = TRUE, strictFilename=FALSE)

Arguments

coordinates | a 2-column numeric matrix of coordinates
df | a data frame object
file | a file name of maximum 8 chararacters, numbers or the underscore symbol to be written, omitting the extensions *.shp, *.shx and *.dbf, which are added in the function
factor2char | logical, default TRUE, convert factor columns to character
strictFilename | if TRUE, impose file basename length limit of 8 characters

Details

The function calls code from shapelib to write a point-type shapefile (both the geometry file with a *.shp extension, the index file with a *.shx extension, and the database file with a *.dbf extension - see write.dbf for details of the data frame export within this function.

Value

no return value.

Note

From maptools 0.4-7, this function is placed in the user-visible namespace on a trial basis, and reports of any malfunction should be sent to the package maintainer, Roger Bivand Roger.Bivand@nhh.no. It is likely that this function and its arguments will be changed.

Author(s)

Nicholas J. Lewin-Koh, modified by Roger Bivand; shapelib by Frank Warmerdam

References

http://shapelib.maptools.org/

See Also

write.dbf

Examples

balt_orig <- readShapePoints(system.file("shapes/baltim.shp", package="maptools")[1])
plot(balt_orig)
balt_cheap <- balt_orig[balt_orig$PRICE < 40,]
file <- tempfile("")
write.pointShape(coordinates=coordinates(balt_cheap), df=as(balt_cheap, "data.frame"), file)
getinfo.shape(paste(file, ".shp", sep=""))
balt_new <- readShapePoints(paste(file, ".shp", sep=""))
plot(balt_new, col="red", pch=16, add=TRUE)


Write a polygon-type shapefile

Description

The function calls code from shapelib to write a polygon-type shapefile from an S3 polylist object.

Usage

write.polylistShape(polylist, df, file, factor2char = TRUE, strictFilename=FALSE, force = FALSE)

Arguments

polylist | list of polygons of class polylist
df | a data frame object
file | a file name of maximum 8 chararacters, numbers or the underscore symbol to be written, omitting the extensions *.shp, *.shx and *.dbf, which are added in the function
factor2char | logical, default TRUE, convert factor columns to character
strictFilename | if TRUE, impose file basename length limit of 8 characters
force | default FALSE, can be set to TRUE to try to force malformed polylist objects to some reasonable form that will not cause both maptools and R to crash. Because polylist objects are old-style rather than new-style classes, it is possible to crash R by trying to write malformed objects. Attempts are made to check for known problems, but using polylist objects rather than sp SpatialPolygons objects is known to be more risky.

Details

The function calls code from shapelib to write a polygon-type shapefile (both the geometry file with a *.shp extension, the index file with a *.shx extension, and the database file with a *.dbf extension - see write.dbf for details of the data frame export within this function.

Value

no return value.

Note

From maptools 0.4-7, this function is placed in the user-visible namespace on a trial basis, and reports of any malfunction should be sent to the package maintainer, Roger Bivand Roger.Bivand@nhh.no. It is likely that this function and its arguments will be changed.

Author(s)

Nicholas J. Lewin-Koh, modified by Roger Bivand; shapelib by Frank Warmerdam

References

http://shapelib.maptools.org/

See Also

write.pointShape, write.dbf

Examples

col_orig <- read.shape(system.file("shapes/columbus.shp", package="maptools")[1])
mappolys <- Map2poly(col_orig, as.character(col_orig$att.data$NEIGNO))
plot(mappolys)
col_df <- col_orig$att.data
col_cheap <- subset(mappolys, col_df$HOVAL < 34)
col_df_cheap <- subset(col_df, col_df$HOVAL < 34)
file <- tempfile("")
write.polylistShape(col_cheap, col_df_cheap, file)
getinfo.shape(paste(file, ".shp", sep=""))
col_new <- read.shape(paste(file, ".shp", sep=""))
mappolys <- Map2poly(col_new, as.character(col_new$att.data$NEIGNO))
plot(mappolys, border="red", add=TRUE)