degreelen, epsg work, with some reasonably good documentation

DESCRIPTION

@@ -1,11 +1,21 @@
 Package: deform
 Type: Package
-Title: What the package does (short line)
-Version: 1.0
+Title: Crustal deformation tools in R
+Version: 0.0-1
 Date: 2014-10-20
-Author: Who wrote it
-Maintainer: Who to complain to <[email protected]>
-Description: More about what it does (maybe more than one line)
-License: What Licence is it under ?
-Imports: Rcpp (>= 0.11.2)
+Author: Andrew J. Barbour; some of these tools were
+    adapted from Matlab code by Brendan Crowell (UW),
+    and Francois Beauducel (IPGP).
+Maintainer: A.J.Barbour <[email protected]>
+Description: Crustal deformation tools.
+License: GPL3
+Imports:
+    Rcpp (>= 0.11.2)
+Depends:
+    RCurl,
+    rgdal,
+    sp,
+    maptools
+Suggests:
+    maps
 LinkingTo: Rcpp

NAMESPACE

@@ -1,3 +1,10 @@
-useDynLib(deform)
+# Generated by roxygen2 (4.0.1): do not edit by hand
+
+export(D_project)
+export(agency.root)
+export(degreelen)
+export(ellipsoid)
+export(getEPSG)
 exportPattern("^[[:alpha:]]+")
-importFrom(Rcpp, evalCpp)
+importFrom(Rcpp,evalCpp)
+useDynLib(deform)

NOBUILD/Crowell.m/strain.m

@@ -25,11 +25,11 @@
         for k=1:length(lon)   
             dx1 = (lon(k)-LON(i,j))*llon(i,j);
             dx2 = (lat(k)-LAT(i,j))*llat(i,j);
-            d = (dx1^2+dx2^2)^0.5/1000;%distance to each grid point
-            u=[u;e(k)/1000;n(k)/1000];%displacements or velocities
-            G=[G;1 0 dx1 dx2 0 0;0 1 0 0 dx1 dx2];%green's function representation of translation plus velocity gradient tensor components
-            W(2*(k-1)+1,2*(k-1)+1)=exp(-d^2/2/aa^2);%Weighting matrix, for east
-            W(2*(k-1)+2,2*(k-1)+2)=exp(-d^2/2/aa^2);%Weighting matrix, for north
+            d = (dx1^2 + dx2^2)^0.5/1000; %distance to each grid point
+             u =[u;e(k)/1000;n(k)/1000]; %displacements or velocities
+            G = [G;1 0 dx1 dx2 0 0;0 1 0 0 dx1 dx2]; %green's function representation of translation plus velocity gradient tensor components
+            W(2*(k-1)+1,2*(k-1)+1)=exp(-d^2/2/aa^2); %Weighting matrix, for east
+            W(2*(k-1)+2,2*(k-1)+2)=exp(-d^2/2/aa^2); %Weighting matrix, for north
         end
 
         [S] = lsqlin(G'*W*G,G'*W*u,[],[]);%solve for strain

R/deform-package.R

@@ -0,0 +1,9 @@
+#' The \code{deform} package:
+#' crustal deformation tools in R.
+#' 
+#' @name deform
+#' @docType package
+#' @useDynLib deform
+#' @exportPattern "^[[:alpha:]]+"
+#' @importFrom Rcpp evalCpp
+NULL

R/degreelen.R

@@ -1,29 +1,57 @@
-#degreelen.m
-#This takes in the latitude and computes the length of 1 degree of latitude
-#and longitude
-#lat is entered in degrees
-
-# ITRF ellipsoid (GRS80), or WGS84 ??
-degreelen <- function(lat, Ea=6378137.0, Eb=6356752.3142){
+#' Compute the length of 1-degree of latitude and longitude
+#' @export
+#' @param lat numeric; degrees
+#' @param model (optional) character; which ellipsoidal correction should be used?
+#' @author A.J.Barbour, adapted from \code{degreelen.m} by B.W.Crowell
+#' @references \url{http://en.wikipedia.org/wiki/World_Geodetic_System}
+#' @examples
+#' lats <- 30:32
+#' # Lon should be about 111.4 km, and Lat 96 - 94 km
+#' degreelen(lats)
+#' 
+#' # In most cases there is no discerible effect
+#' # caused by varying the semi-major axis.
+#' #   Not much difference according to this machine:
+#' all.equal(degreelen(lats,"IRTF"), degreelen(lats,"WGS84"), tolerance = .Machine$double.eps)
+#' 
+#' #   unless the precision of the input is in micrometers, 
+#' #   which is essentially impossible to achieve:
+#' options(digits=13)
+#' lats <- 30.1111111123
+#' rbind(degreelen(lats,"IRTF"),degreelen(lats,"WGS84"))
+degreelen <- function(lat, model){
   #
-  #http://surveying.wb.psu.edu/sur162/control/surfaces.htm
-  Fe. <- 1 - Eb/Ea # flattening
-  EE. <- 2*Fe. - Fe. * Fe. # eccentricity
+  pic <- pi/180
+  latrad <- as.numeric(lat) * pic
+  #
+  if (missing(model)) model <- "IRTF"
+  Ell <- ellipsoid(model)
+  Ea. <- Ell[["major"]] # major axis
+  Fe. <- Ell[["flat"]]  # flattening
+  EE. <- Ell[["ecc"]]   # eccentricity
   #
-  pic = pi/180
-  latrad = lat * pic
+  #http://surveying.wb.psu.edu/sur162/control/surfaces.htm
   slp <- sin(latrad)
   slpsq <- slp * slp
   #
-  degDist <- 1 - EE. * slpsq
+  degDist <- 1 - EE. * slp * slp
   degDistsq <- sqrt(degDist)
   #
-  N. <- Ea / degDistsq
-  M. <- Ea * cos(latrad) * (1 - EE.) / degDistsq / degDist ## equiv to ^1.5
+  N. <- Ea. / degDistsq
+  M. <- Ea. * cos(latrad) * (1 - EE.) / degDistsq / degDist ## equiv to ^1.5
   #
-  toret <- zapsmall(cbind(N., M.) * pic)
-  rownames(toret) <- lat
-  return(toret)
+  return(cbind(Lat=lat, (cbind(Meters.lon=N., Meters.lat=M.) * pic)))
+}
+
+#' @rdname degreelen
+#' @export
+ellipsoid <- function(model=c("IRTF","WGS84")){
+  # ITRF ellipsoid (GRS80) is default
+  mdl <- match.arg(model)
+  message(mdl)
+  Ea <- 6378137.0
+  Eb <- 6356752.0 + switch(mdl, IRTF=0.314140, WGS84=0.314245)
+  Flat <- 1 - Eb/Ea       # flattening
+  Ecc <- Flat*(2 - Flat)  # eccentricity
+  return(data.frame(model=mdl, major=Ea, minor=Eb, flat=Flat, ecc=Ecc))
 }
-##
-##

R/epsg.R

@@ -0,0 +1,99 @@
+#' Functions to work with EPSG projection specifications
+#' @name epsg
+#' @aliases EPSG
+#' @param agency character;
+#' @param ref.id integer;
+#' @param verbose logical;
+#' @param Df an object to be coerced into a \code{\link{data.frame}}
+#' @param coords character; the names of the coordinates in \code{Df}
+#' @param p4,p4old; these are the proj4 strings for the new (\code{p4})
+#'  and old (\code{p4old}) data.  More specifically, \code{p4} will be
+#'  the new projection (set by \code{\link{CRS}}), and  
+#'  \code{p4old} is the old projection (set by \code{\link{proj4string}}).
+#' @param ... additional parameters to \code{\link{elide}}
+#' @references 
+#' A list of EPSG projections:
+#'   \url{http://spatialreference.org/ref/epsg/}
+#' @references 
+#' A list of PROJ commands:
+#'   \url{https://trac.osgeo.org/proj/wiki/GenParms}
+#' @examples
+#' \dontrun{
+#' library(maps)
+#' library(maptools)
+#' library(datasets)
+#' data(state)
+#' states <- data.frame(state.x77, state.center)
+#' # HI/AK are in the middle of the ocean, so take them out:
+#' states <- states[states$x > -121,]
+#' 
+#' # uses getEPSG to project if p4/p4old are missing
+#' d1 <- D_project(states, c("x","y"))
+#' 
+#' # set p4
+#' p <- getEPSG()
+#' d2 <- D_project(states, c("x","y"), p)
+#' all.equal(d1,d2) # TRUE
+#' 
+#' # reproject into UTM-11 in kilometers
+#' d2u <- spTransform(d2, CRS("+proj=utm +zone=11 +ellps=WGS84 +units=km"))
+#' 
+#' # change the projection
+#' pold <- p
+#' p <- getEPSG(2153)  # UTM-11 in meters
+#' d3 <- D_project(states, c("x","y"), p, pold)
+#' 
+#' layout(matrix(1:6,ncol=3))
+#' us <- map("usa",plot=FALSE)
+#' plot(d1, axes=TRUE); title("d1"); lines(us)
+#' plot(d2, axes=TRUE);title("d2 (== d1)"); lines(us)
+#' plot(d2u, axes=TRUE);title("d2 - utm11/km")
+#' plot(d3, axes=TRUE);title("d3 - utm11/m")
+#' #
+#' # rotate with maptools (can also do this using ...)
+#' plot(elide(d2u, rotate=45));title("d2 at 45") # clockwise
+#' plot(elide(d2u, rotate=-45));title("d2 at -45") # counter-clockwise
+#' #
+#' }
+NULL
+
+#' @rdname epsg
+#' @export
+agency.root <- function(agency=NULL){
+  agency <- match.arg(agency,c("epsg","esri","iau2000","sr-org"))
+  return(sprintf("http://spatialreference.org/ref/%s",agency))
+}
+
+#' @rdname epsg
+#' @export
+getEPSG <- function(ref.id=4326, verbose=TRUE, agency=NULL){
+  ref.id <- as.integer(ref.id)
+  #http://spatialreference.org/ref/epsg/2805/proj4/
+  #http://spatialreference.org/ref/epsg/4326/proj4/
+  ref.url <- sprintf("%s/%s/proj4/", agency.root(agency), ref.id)
+  if (verbose) message(ref.url)
+  stopifnot(RCurl::url.exists(ref.url))
+  prj <- RCurl::getURL(ref.url)
+  return(prj)
+}
+
+#' @rdname epsg
+#' @export
+D_project <- function(Df, coords=names(Df)[1:2], p4, p4old, verbose=TRUE, ...){
+  Df <- data.frame(Df)
+  stopifnot(length(coords)>=2)
+  #
+  # proj4 string
+  if (missing(p4)) p4 <- getEPSG(verbose=verbose)
+  if (missing(p4old)) p4old <- p4
+  eprj <- CRS(p4)
+  stopifnot(exists("eprj"))
+  ##
+  coordinates(Df) <- coords
+  proj4string(Df) <- p4old
+  #
+  toret <- try(elide(spTransform(Df, CRSobj=eprj), ...))
+  stopifnot(!inherits(toret,"try-error"))
+  return(toret)
+}
+##

man/deform.Rd

@@ -0,0 +1,12 @@
+% Generated by roxygen2 (4.0.1): do not edit by hand
+\docType{package}
+\name{deform}
+\alias{deform}
+\alias{deform-package}
+\title{The \code{deform} package:
+crustal deformation tools in R.}
+\description{
+The \code{deform} package:
+crustal deformation tools in R.
+}
+

man/degreelen.Rd

@@ -0,0 +1,41 @@
+% Generated by roxygen2 (4.0.1): do not edit by hand
+\name{degreelen}
+\alias{degreelen}
+\alias{ellipsoid}
+\title{Compute the length of 1-degree of latitude and longitude}
+\usage{
+degreelen(lat, model)
+
+ellipsoid(model = c("IRTF", "WGS84"))
+}
+\arguments{
+\item{lat}{numeric; degrees}
+
+\item{model}{(optional) character; which ellipsoidal correction should be used?}
+}
+\description{
+Compute the length of 1-degree of latitude and longitude
+}
+\examples{
+lats <- 30:32
+# Lon should be about 111.4 km, and Lat 96 - 94 km
+degreelen(lats)
+
+# In most cases there is no discerible effect
+# caused by varying the semi-major axis.
+#   Not much difference according to this machine:
+all.equal(degreelen(lats,"IRTF"), degreelen(lats,"WGS84"), tolerance = .Machine$double.eps)
+
+#   unless the precision of the input is in micrometers,
+#   which is essentially impossible to achieve:
+options(digits=13)
+lats <- 30.1111111123
+rbind(degreelen(lats,"IRTF"),degreelen(lats,"WGS84"))
+}
+\author{
+A.J.Barbour, adapted from \code{degreelen.m} by B.W.Crowell
+}
+\references{
+\url{http://en.wikipedia.org/wiki/World_Geodetic_System}
+}
+