R/orthproj.R

In geo: Draw and Annotate Maps, Especially Charts of the North Atlantic

orthproj <-
function(pts, curve)
{
        pts1 <- lambert(pts$lat, pts$lon, 65, -18, 65)
        curve1 <- lambert(curve$lat, curve$lon, 65, -18, 65)
        pts$x <- pts1$x/1.852
        pts$y <- pts1$y/1.852
        curve$x <- curve1$x/1.852
        curve$y <- curve1$y/1.852
        pardist <- perdist <- rep(0, length(pts$lat))
        x <- .C("Curvedist", PACKAGE = "geo", 
                as.double(curve$x),
                as.double(curve$y),
                as.double(curve$dist),
                as.integer(length(curve$x)),
                as.double(pts$x),
                as.double(pts$y),
                as.double(perdist),
                as.double(pardist),
                as.integer(length(pts$y)))
        pardist <- x[[7]]
        perdist <- x[[8]]
        # Points inside curve get negative perdist.
        i <- geoinside(pts, reg = curve, option = 0, robust = F)
        perdist[i] <-  - perdist[i]
        return(list(pardist = pardist, perdist = perdist))
}