[R-sig-Geo] Maps across the dateline

[Roger Bivand]

On Thu, 19 Jan 2006 White.Denis at epamail.epa.gov wrote:

While there are plenty of good answers to this question, the most direct 
one - given that the maps package is maintained from New Zealand, is to 
use the world2 database it provides:

?world2

which runs 0-360 (roughly), and should be good enough resolution for a 
whole-Pacific map.

Roger

> 
> 
> >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)
> 
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> 

-- 
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no