[R-sig-Geo] Maps across the dateline

[White.Denis at epamail.epa.gov]

>I'm trying to plot some
>data for the Pacific Ocean and want to get a map onto my plot as
>well.
>I'm using the function 'map()' from library(maps). However, I cannot
>get it to plot anything across the dateline since it requires that
>longitudes between 180 and 360 degrees east to be written as negative
>degrees east of Greenwich (e.g. -180 to 0 degrees east). The
>longitudes
>I want to plot are 135 to 255 degrees east but the function map()
>wants
>that to be written 135 to 180 degrees east and -180 to -105 degrees
>east, so I would have to divide the grid up in two. But even if I do
>that I can't plot it together because -180 doesn't come after 180.
>
>Does anyone have a way to plot something across the dateline? I only
>want a map of the Pacific, not the rest of the world.
>
>Thanks,
>Anne Hertel

If you use an appropriate map projection centered at 180/-180, for
example, then you should be able to plot your data and islands or
latitude/longitude lines as desired.  See the following example for
plotting some random points and the latitude and longitude graticule
lines using one particular map projection.  This projection, the Lambert
azimuthal, is an equal area projection but may not be the most
appropriate for your needs.  There are web sites that can help to choose
a projection if you need it.

# plot data across 180/-180

# Lambert azimuthal equal area map projection
lambaz <- function (lon, lat, clon=NULL, clat=NULL)
  # Arguments are vectors of longitude and latitude,
  # that is, geographic coordinates in decimal degrees.
  # Projection center is (clon, clat), which, if null, is
  # set to the midpoint of the input in lat-lon space.
  # Line identifiers in lon[is.na (lat)] retained in x.
  # Returns a two column matrix with column names "x" and "y",
  # in units of kilometers.  This is Lambert's azimuthal
  # equal area map projection using a spherical earth model.
  # (ref JP Snyder, USGS Prof Paper 1395, p 182 ff).
{
    R <- 6371 # authalic radius of Clarke 1866 rounded to km
    if (is.null (clon) || is.null (clat))
    {
        rlat <- range (lat, na.rm=TRUE)
        rlon <- range (lon, na.rm=TRUE)
        clat <- diff (rlat) / 2 + rlat[1]
        clon <- diff (rlon) / 2 + rlon[1]
    }
    dlon <- (lon - clon) * pi / 180
    cosdlon <- cos (dlon)
    sinlat <- sin (lat * pi / 180)
    coslat <- cos (lat * pi / 180)
    sinclat <- sin (clat * pi / 180)
    cosclat <- cos (clat * pi / 180)
    kp <- 1 + sinclat*sinlat + cosclat*coslat*cosdlon
    kp <- sqrt (2 / kp)
    x <- R * kp * coslat * sin (dlon)
    y <- R * kp * (cosclat*sinlat - sinclat*coslat*cosdlon)
    x[is.na (lat)] <- lon[is.na (lat)]
    xy <- cbind (x=x, y=y)
    xy
}

# generate latitude and longitude graticules
graticule <- function (minlon=-180, maxlon=180,
    minlat=-90, maxlat=90, ginc=10, dinc=2)
{
    dlon <- maxlon - minlon
    dlat <- maxlat - minlat
    nr <- (dlon / dinc + 1) * (dlat / ginc + 1) +
        (dlat / dinc + 1) * (dlon / ginc + 1) +
        dlon / ginc + dlat / ginc + 2
    geo <- matrix (0, nrow=nr, ncol=2)
    i <- n <- 0
    for (lon in seq (minlon, maxlon, by=ginc))
    {
        i <- i + 1
        n <- n + 1
        geo[i, 1] <- n
        geo[i, 2] <- NA
        for (lat in seq (minlat, maxlat, by=dinc))
        {
            i <- i + 1
            geo[i, 1] <- lon
            geo[i, 2] <- lat
        }
    }
    for (lat in seq (minlat, maxlat, by=ginc))
    {
        i <- i + 1
        n <- n + 1
        geo[i, 1] <- n
        geo[i, 2] <- NA
        for (lon in seq (minlon, maxlon, by=dinc))
        {
            i <- i + 1
            geo[i, 1] <- lon
            geo[i, 2] <- lat
        }
    }
    geo
}

# a graticule for central pacific
east.part <- graticule (minlon=-180, maxlon=-150,
                        minlat=-30, maxlat=30)
west.part <- graticule (minlon= 150, maxlon= 180,
                        minlat=-30, maxlat=30)
grat.geo <- rbind (east.part, west.part)
grat.lam <- lambaz (grat.geo[,1], grat.geo[,2], clon=180, clat=0)

# generate random points in central pacific
east.pts <- cbind (runif (100, min=-180, max=-150),
                   runif (100, min=-30, max=30))
west.pts <- cbind (runif (100, min= 150, max= 180),
                   runif (100, min=-30, max=30))
pts.geo <- rbind (east.pts, west.pts)
pts.lam <- lambaz (pts.geo[,1], pts.geo[,2], clon=180, clat=0)

# plot in geographic coordinate space (lat, lon)
plot.new ()
plot.window (range (grat.geo[,1], na.rm=TRUE),
             range (grat.geo[,2], na.rm=TRUE), asp=1)
lines (grat.geo)
points (pts.geo, pch=20)

# plot in Lambert space
plot.new ()
plot.window (range (grat.lam[,1], na.rm=TRUE),
             range (grat.lam[,2], na.rm=TRUE), asp=1)
lines (grat.lam)
points (pts.lam, pch=20)