Re: [R-sig-Geo] Good projection for N/S America?

[Roger Bivand]

On Mon, 9 Apr 2007 [EMAIL PROTECTED] wrote:

> For preview graphics and for large areas such as continents, large
> countries, hemispheres, or the whole earth, spherical projections are
> often adequate.  I can provide some of the ones I have used.  For
> detailed work at sites and small areas, ellipsoidal projections such as
> UTM are usually used, and then the coding gets more complicated with
> choices of datums and so forth.
> 

The attached script shows how to do the interrupted sinusoidal projection 
using spTransform in rgdal, for the whemi.lin data posted with the 
free-standing functions by Denis White a couple of days ago. Once the 
lines are converted into SpatialLines objects, the rest is robust and 
simple, as is the use of gridlines() in sp. The one catch is calculating 
the offset, here in an x_0= offset along the Equator in metres between the 
two central longitude values. The output is attached as a PNG image. The 
point about the sp objects is that they contain enough metadata (here a 
PROJ.4 projection description) to let them be moved to other R packages or 
external software. 

The half-dozen basic projections are easy to specify in PROJ.4, for 
example from the geotiff list:

http://www.remotesensing.org/geotiff/proj_list/

which is what I used here. The other projections mentioned are:

Lambert Cylindrical Equal Area "+proj=cea +lon_0=-80"

Lambert Azimuthal Equal Area "+proj=laea +lat_0=0 +lon_0=-80"

while the Northern hemisphere sinusoidal is:

Sinusoidal "+proj=sinu +lon_0=-100"

So I'd argue that PROJ.4 projection descriptions are not difficult to use, 
and with sp objects, do stay stuck to the data (has anyone else ever 
forgotten what projection was used when revisiting data, not just me?).

Using the maptools map2SpatialLines() interface function, or the Rgshhs() 
interface to GSHHS shorelines, even getting the lines is quite easy, 
qualified by clipping and bounding box issues in extremities for 
projection from geographical coordinates.

Of course, it would help to have MacOS X and selected Linux binaries of 
rgdal, we're very lucky that Uwe Ligges is so helpful with the Windows 
binaries.

Roger

> [EMAIL PROTECTED] wrote on 2007-04-08 07:56:03:
> 
> > Denis,
> >
> > That's really useful. It occurs to me that we only really need a
> > half-dozen basic projections to cover 90% of user cases. Perhaps these
> > could be incorporated into the 'sp' group somewhere and relieve the
> > dependence on proj4. (It could be packaged separately for R  for the
> > other 10% of cases where its needed.)
> >
> > THK
> >
> > On 4/6/07, [EMAIL PROTECTED] <[EMAIL PROTECTED]>
> wrote:
> > > Thanks, Roger.  There was a request to see the R code for these
> figures.
> > > Attached is the script for the second PDF file plus the input
> boundary
> > > file I used for the hemisphere.  The three projection functions are
> for
> > > simple spherical, rather than ellipsoidal, models of the earth.  The
> > > graticule generating function could be more elegant.  I'm not yet up
> to
> > > speed with sp and the many new spatial capabilities in R so please
> > > excuse the old style "lines()" format encoding and graphics.
> > >
> > > Tim, I don't know whether proj4 could do the interrupted sinusoidal.
> > >
> > > (See attached file: whemi.projs.r)(See attached file: whemi.lin)
> > >
> > > [EMAIL PROTECTED] wrote on 2007-04-06 04:51:53:
> > >
> > > > Since this topic is of general interest, I've made an exception
> and
> > > > allowed (this once!) a posting of more than 200K. In general, if
> > > graphics
> > > > are big, please consider either an alternative device (png is
> often
> > > OK),
> > > > or posting just a URL to the real file.
> > > >
> > > > With apologies to list members on dial-up connections in the
> field,
> > > >
> > > > Roger
> > > >
> > > >
> > > > On Thu, 5 Apr 2007 [EMAIL PROTECTED] wrote:
> > > >
> > > > > Yes, for many uses that is my choice also.  For the conterminous
> US
> > > for
> > > > > example, the Lambert azimuthal has lower mean distortion than
> the
> > > > > commonly used standard projection, the Albers conical equal
> area,
> > > > > although Albers was chosen by USGS as a standard because of
> lower
> > > > > extreme distortion than many other possible projections.
> > > > >
> > > > > For our hemispherical application, because we were gridding the
> > > data, we
> > > > > wanted parallels of latitude to be parallel in the projected
> > > coordinate
> > > > > space, which we wouldn't get with the Lambert azimuthal.
> > > > >
> > > > > (See attached file: whemi.projs.pdf)
> > > > >
> > > > > Tim Keitt <[EMAIL PROTECTED]> wrote on 2007-04-05 10:56:09:
> > > > >
> > > > > > Thanks. My application is not that demanding. Really, I just
> want
> > > it
> > > > > > to look reasonable. My plan is to lay out the postings in the
> > > > > > projected coordinates and then back transform into geographic
> > > > > > coordinates for analysis. I tried lots of projections and
> found
> > > > > > Lamberts Azimuthal Equal Area to be quite good. I like the
> look of
> > > the
> > > > > > Azimuthal Equidistant better, but figured equal area was a
> good
> > > > > > choice.
> > > > > >
> > > > > > THK
> > > > > >
> > > > > > On 4/4/07, [EMAIL PROTECTED]
> > > <[EMAIL PROTECTED]>
> > > > > wrote:
> > > > > > > Tim,
> > > > > > >
> > > > > > > It depends on which kind of distortion is of most concern.
> For
> > > many
> > > > > > > types of extensive data, especially counts, for example, the
> > > equal
> > > > > area
> > > > > > > property is desirable.  We used the Lambert cylindrical
> equal
> > > area
> > > > > > > projection with standard parallels of +/- 30 degrees for
> some
> > > > > western
> > > > > > > hemispherical work, see reference below.  (The center
> longitude
> > > > > could be
> > > > > > > -80 west, but that is less important than the choice of
> > > parallels.)
> > > > > > >
> > > > > > > Before falling back on the Lambert as an easy to use
> projection,
> > > I
> > > > > tried
> > > > > > > to get several ESRI products to implement an interrupted
> > > projection
> > > > > > > using the sinusoidal projection, in part for reasons given
> in
> > > the
> > > > > second
> > > > > > > reference.  I used a separate center longitude for north and
> > > south
> > > > > of
> > > > > > > the equator and the appearance is certainly more
> satisfactory
> > > than
> > > > > the
> > > > > > > Lambert in my opinion.  I'll attach a PDF of an illustration
> of
> > > this
> > > > > > > approach generated in R that I hope you will get but not the
> > > rest of
> > > > > the
> > > > > > > list unfortunately.  I can send PDFs of the references also
> if
> > > > > needed.
> > > > > > >
> > > > > > > Denis
> > > > > > >
> > > > > > > Lawler JJ, White D, Neilson RP, Blaustein AR.  2006.
> Predicting
> > > > > > > climate-induced range shifts: model differences and model
> > > > > reliability.
> > > > > > > Global Change Biology 12:1568-1584.
> > > > > > >
> > > > > > > White D.  2006.  Display of pixel loss and replication in
> > > > > reprojecting
> > > > > > > raster data from the sinusoidal projection.  Geocarto
> > > International
> > > > > > > 21(2):19-22.
> > > > > > >
> > > > > > > (See attached file: whemi.sinus.pdf)
> > > > > > >
> > > > > > > [EMAIL PROTECTED] wrote on 2007-04-04
> > > 12:17:39:
> > > > > > >
> > > > > > > > Anyone know of a particularly good map projection for
> showing
> > > all
> > > > > of
> > > > > > > > North and South America without too much distortion?
> > > > > > > >
> > > > > > > > THK
> > > > > > > >
> > > > > > > > --
> > > > > > > > Timothy H. Keitt, University of Texas at Austin
> > > > > > > > Contact info and schedule at
> http://www.keittlab.org/tkeitt/
> > > > > > > > Reprints at http://www.keittlab.org/tkeitt/papers/
> > > > > > > > ODF attachment? See http://www.openoffice.org/
> > > > > > > >
> > > > > > > > _______________________________________________
> > > > > > > > R-sig-Geo mailing list
> > > > > > > > R-sig-Geo@stat.math.ethz.ch
> > > > > > > > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> > > > > > >
> > > > > >
> > > > > >
> > > > > > --
> > > > > > Timothy H. Keitt, University of Texas at Austin
> > > > > > Contact info and schedule at http://www.keittlab.org/tkeitt/
> > > > > > Reprints at http://www.keittlab.org/tkeitt/papers/
> > > > > > ODF attachment? See http://www.openoffice.org/
> > > >
> > > > --
> > > > 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: [EMAIL PROTECTED]
> > > >
> > > > _______________________________________________
> > > > R-sig-Geo mailing list
> > > > R-sig-Geo@stat.math.ethz.ch
> > > > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> > >
> >
> >
> > --
> > Timothy H. Keitt, University of Texas at Austin
> > Contact info and schedule at http://www.keittlab.org/tkeitt/
> > Reprints at http://www.keittlab.org/tkeitt/papers/
> > ODF attachment? See http://www.openoffice.org/
> >
> > _______________________________________________
> > R-sig-Geo mailing list
> > R-sig-Geo@stat.math.ethz.ch
> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> 
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo@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: [EMAIL PROTECTED]


whemi <- read.table ("whemi.lin", header=TRUE)
# make the input object into two lists of lines split at Equator
whemiNA <- is.na(whemi$x) | is.na(whemi$y)
whemiCS <- cumsum(whemiNA)
whemi1 <- cbind(whemi, whemiCS)
whemi2 <- whemi1[!whemiNA,]
whemiL <- split(whemi2, whemi2$whemiCS)
whemiL_NS <- sapply(whemiL, function(x) { r <- range(x$y) ;
    (r[1] >= 0 && r[2] <= 0) || (r[1] <= 0 && r[2] >=0) })
whemiL1 <- whemiL[!whemiL_NS]
whemiL2 <- whemiL[whemiL_NS]
whemiL2a <- vector(mode="list", length=length(whemiL2)*2)
ii <- 1
for (i in seq(along=whemiL2)) {
    cy <- which(whemiL2[[i]]$y == 0.0)
    n <- length(whemiL2[[i]]$y)
    whemiL2a[[ii]] <- whemiL2[[i]][1:cy,]
    ii <- ii+1
    whemiL2a[[ii]] <- whemiL2[[i]][cy:n,]
    ii <- ii+1
}
whemiL3 <- c(whemiL1, whemiL2a)
whemiNorth <- sapply(whemiL3, function(x) any(x$y > 0))
whemiN <- whemiL3[whemiNorth]
whemiS <- whemiL3[!whemiNorth]
# convert N and S line lists into SpatialLines objects
library(sp)
whemiNLs <- lapply(whemiN, function(x) Lines(list(Line(cbind(x$x, x$y))),
    ID=as.character(x$whemiCS[1])))
whemiSLs <- lapply(whemiS, function(x) Lines(list(Line(cbind(x$x, x$y))),
    ID=as.character(x$whemiCS[1])))
whemiN_SL <- SpatialLines(whemiNLs,
    proj4string=CRS("+proj=longlat +datum=WGS84"))
whemiS_SL <- SpatialLines(whemiSLs,
    proj4string=CRS("+proj=longlat +datum=WGS84"))
# add grids as SpatialLines objects
grdN <- gridlines(whemiN_SL, easts=seq(-170,-30,10), norths=seq(0,80,10),
    ndiscr=500)
grdS <- gridlines(whemiS_SL, easts=seq(-170,-30,10), norths=seq(-50,0,10),
    ndiscr=500)
library(rgdal)
# compute offset for interrupted sinusoidal projection
onept <- matrix(c(-80,0), nrow=1)
oneSP <- SpatialPoints(onept, proj4string=CRS("+proj=longlat +datum=WGS84"))
oneSPN <- spTransform(oneSP, CRS("+proj=sinu +lon_0=-100"))
oneSPS <- spTransform(oneSP, CRS("+proj=sinu +lon_0=-60"))
offset <- coordinates(oneSPN)[,1] - coordinates(oneSPS)[,1]
# project south SpatialLines objects with offset
grdS.sinus <- spTransform(grdS, CRS(paste("+proj=sinu +lon_0=-60 +x_0=",
    offset, sep="")))
whemiS.sinus <- spTransform(whemiS_SL,
    CRS(paste("+proj=sinu +lon_0=-60 +x_0=", offset, sep="")))
# and north SpatialLines objects without
grdN.sinus <- spTransform(grdN, CRS("+proj=sinu +lon_0=-100"))
whemiN.sinus <- spTransform(whemiN_SL, CRS("+proj=sinu +lon_0=-100"))
# assemble the joint bounding box
bb_all <- t(apply(cbind(bbox(whemiS.sinus), bbox(whemiN.sinus)), 1, range))
plot(whemiN.sinus, xlim=bb_all[1,], ylim=bb_all[2,])
plot(whemiS.sinus, add=TRUE)
plot(grdS.sinus, add=TRUE, col="red", lty=2)
plot(grdN.sinus, add=TRUE, col="red", lty=2)

whemi.png
Description: PNG image

_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo