Introduction

Through the early months of 2015, efforts were made to incorporate PROJ.4 projections into the oce R package. An earlier attempt had relied on the proj4 R package, but this package has the disadvantage of reporting an error if any point in a polygon cannot be projected, so it had to be called pointwise on coastline vectors, at significant time cost (minutes for many practical cases). The next step was to incorporate PROJ.4 code into oce itself. This was used in oce version 0.9-17. However, a CRAN-core member suggested removing the PROJ.4 code from oce when 0.9-17 was being vetted, and so later versions of ‘oce’ rely on the ‘rgdal’ package.

Oce only provides those PROJ.4 projections that have inverses. This is because oce needs to do inverse projections for some of its plotting operations. Also note that a handful of projections have not been incorporated, because they cause errors of various sorts; see oce issue 635 for more on this.

I will not explain the methodology of projection here, nor the details of any given projection. Readers should consult the help for the oce function mapPlot() and the many websites detailing software that derives from PROJ.4. In due course, they will be able to consult my book.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
library(oce)
data(coastlineWorld)
cl45 <- coastlineCut(coastlineWorld, lon_0=-45)

par(mar=rep(2, 4))
line <- 0.25
pcol <- "blue"
ecol <- "red"
font <- 2
col <- "lightgray"

p <- "+proj=aea +lat_1=10 +lat_2=60 +lon_0=-40"
mapPlot(coastlineWorld, projection=p, longitudelim=c(-80,0), latitudelim=c(0,45), col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
4
p <- "+proj=aeqd +lon_0=-45"
mapPlot(coastlineWorld, projection=p,
        longitudelim=c(-80,0), latitudelim=c(0,60), col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
4
5
p <- "+proj=aitoff +lon_0=-45"
mapPlot(coastlineWorld, projection=p,
        longitudelim=c(-80,0), latitudelim=c(0,60), col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)
mtext("Problem at top, unless coastlineCut() used", line=line, adj=0, col=ecol, font=font)

center

1
2
3
4
5
p <- "+proj=aitoff +lon_0=-45"
mapPlot(cl45, projection=p,
        longitudelim=c(-80,0), latitudelim=c(0,60), col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)
mtext("Using coastlineCut()", line=line+0.9, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=bipc"
mapPlot(coastlineWorld, projection=p, longitudelim=c(-80,0), latitudelim=c(70, 110), col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=bonne +lat_1=45"
mapPlot(coastlineWorld, projection=p, col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=cass +lon_0=-45"
mapPlot(coastlineWorld, projection=p, longitudelim=c(-80,0), latitudelim=c(0,60), col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
4
p <- "+proj=cass +lon_0=-45"
mapPlot(cl45, projection=p, longitudelim=c(-80,0), latitudelim=c(0,60), col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)
mtext("Using coastlineCut()", line=line+0.9, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=cc"
mapPlot(coastlineWorld, projection=p, longitudelim=c(-80,0), latitudelim=c(-40,40), col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=cea"
mapPlot(coastlineWorld, projection=p, col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=collg"
mapPlot(coastlineWorld, projection=p, col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=crast"
mapPlot(coastlineWorld, projection=p, col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=eck1"
mapPlot(coastlineWorld, projection=p, col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=eck2"
mapPlot(coastlineWorld, projection=p, col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=eck3"
mapPlot(coastlineWorld, projection=p, col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=eck4"
mapPlot(coastlineWorld, projection=p, col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=eck5"
mapPlot(coastlineWorld, projection=p, col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=eck6"
mapPlot(coastlineWorld, projection=p, col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=eqc"
mapPlot(coastlineWorld, projection=p, col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=euler +lat_1=45 +lat_2=50 +lon_0=-40"
mapPlot(coastlineWorld, projection=p, longitudelim=c(-80, 0), latitudelim=c(0, 60), col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=etmerc +ellps=WGS84 +lon_0=-40"
mapPlot(coastlineWorld, projection=p, longitudelim=c(-80, 0), latitudelim=c(0, 60), col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
4
p <- "+proj=etmerc +ellps=WGS84 +lon_0=-40"
mapPlot(cl45, projection=p, longitudelim=c(-80, 0), latitudelim=c(0, 60), col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)
mtext("Using coastlineCut()", line=line+0.9, adj=1, col=pcol, font=font)

center

1
2
3
4
p <- "+proj=fahey"
mapPlot(coastlineWorld, projection=p, col=col)
mtext("+proj=fahey", line=line, adj=1, col=pcol, font=font)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=fouc"
mapPlot(coastlineWorld, projection=p, col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=fouc_s"
mapPlot(coastlineWorld, projection=p, col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=gall"
mapPlot(coastlineWorld, projection=p, col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=geos +h=1e8"
mapPlot(coastlineWorld, projection=p, col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=gn_sinu +n=6 +m=3"
mapPlot(coastlineWorld, projection=p, col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=gnom +lon_0=-40"
mapPlot(coastlineWorld, projection=p, longitudelim=c(-80,0), latitudelim=c(-30,30), col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=goode"
mapPlot(coastlineWorld, projection=p, col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=hatano"
mapPlot(coastlineWorld, projection=p, col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
4
p <- "+proj=healpix +a=1"
mapPlot(coastlineWorld, projection=p, col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)
mtext("Problem in Canada", line=line, adj=0, col=ecol, font=font)

center

1
2
3
4
p <- "+proj=rhealpix +south_square=0 +north_square=1"
mapPlot(coastlineWorld, projection=p)
mtext(p, line=line, adj=1, col=pcol, font=font)
mtext("Is this an acid trip?", line=line, adj=0, col=ecol, font=font)

center

1
2
3
4
p <- "+proj=igh"
mapPlot(coastlineWorld, projection=p, col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)
mtext("Connection through cutout region", line=line, adj=0, col=ecol, font=font)

center

1
2
3
p <- "+proj=kav5"
mapPlot(coastlineWorld, projection=p, col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=kav7"
mapPlot(coastlineWorld, projection=p, col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=laea +lon_0=-40"
mapPlot(coastlineWorld, projection=p, col=col)
mtext("+proj=laea +lon_0=-40", line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=lonlat"
mapPlot(coastlineWorld, projection=p, col=col)
mtext("+proj=lonlat", line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=latlon"
mapPlot(coastlineWorld, projection=p, col=col)
mtext("+proj=lonlat", line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=lcc +lat_1=30 +lat_2=70 +lon_0=-40"
mapPlot(coastlineWorld, projection=p, longitudelim=c(-80, 0), latitudelim=c(0, 60), col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=leac +lon_0=-40"
mapPlot(coastlineWorld, projection=p, longitudelim=c(-80, 0), latitudelim=c(0, 60), col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=loxim"
mapPlot(coastlineWorld, projection=p, col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=mbt_s"
mapPlot(coastlineWorld, projection=p, col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=mbt_fps"
mapPlot(coastlineWorld, projection=p, col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=mbtfpp"
mapPlot(coastlineWorld, projection=p, col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=mbtfpq"
mapPlot(coastlineWorld, projection=p, col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=mbtfps"
mapPlot(coastlineWorld, projection=p, col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=merc"
mapPlot(coastlineWorld, projection=p, longitudelim=c(-80,0), latitudelim=c(0,45), col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=mil_os"
mapPlot(coastlineWorld, projection=p, longitudelim=c(-80,0), latitudelim=c(70, 120), col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=mill"
mapPlot(coastlineWorld, projection=p, col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=moll"
mapPlot(coastlineWorld, projection=p, col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=murd1 +lat_1=30 +lat_2=60 +lon_0=-40"
mapPlot(coastlineWorld, projection=p, longitudelim=c(-80,0), latitudelim=c(0,45), col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=murd2 +lat_1=30 +lat_2=60 +lon_0=-40"
mapPlot(coastlineWorld, projection=p, longitudelim=c(-80,0), latitudelim=c(0,45), col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=murd3 +lat_1=30 +lat_2=60 +lon_0=-40"
mapPlot(coastlineWorld, projection=p, longitudelim=c(-80,0), latitudelim=c(0,45), col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=natearth"
mapPlot(coastlineWorld, projection=p, col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=nell"
mapPlot(coastlineWorld, projection=p, col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=nell_h"
mapPlot(coastlineWorld, projection=p, col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
4
p <- "+proj=nsper +h=90000000"
mapPlot(coastlineWorld, projection=p, col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)
mtext("Problem in Antarctica", line=line, adj=0, col=ecol, font=font)

center

1
2
3
4
5
6
p <- "+proj=ob_tran"
## mapPlot(coastlineWorld, projection=p, longitudelim=c(-180,-120), latitudelim=c(-50,-20), col=col)
plot(0:1, 0:1, axes=FALSE, type='n', xlab="", ylab="")
box()
mtext(p, line=line, adj=1, col=pcol, font=font)
mtext("Usage?", line=line, adj=0, col=ecol, font=font)

center

1
2
3
4
p <- "+proj=ocea"
mapPlot(coastlineWorld, projection=p, col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)
mtext("Ugly Horizontal Lines", line=line, adj=0, col=ecol, font=font)

center

1
2
3
p <- "+proj=omerc +lat_1=30 +lon_1=-40 +lat_2=60"
mapPlot(coastlineWorld, projection=p, longitudelim=c(-80,0), latitudelim=c(0,60), col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
4
p <- "+proj=ortho"
mapPlot(coastlineWorld, projection=p, col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)
mtext("Problem in Antarctica", line=line, adj=0, col=ecol, font=font)

center

1
2
3
p <- "+proj=pconic +lat_1=20 +lat_2=60 +lon_0=-40"
mapPlot(coastlineWorld, projection=p, longitudelim=c(-80,0), latitudelim=c(0,60), col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=poly +lon_0=-40"
mapPlot(coastlineWorld, projection=p, longitudelim=c(-80,0), latitudelim=c(0,60), col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=putp1"
mapPlot(coastlineWorld, projection=p, col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=putp2"
mapPlot(coastlineWorld, projection=p, col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=putp3"
mapPlot(coastlineWorld, projection=p, col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=putp5"
mapPlot(coastlineWorld, projection=p, col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=putp6"
mapPlot(coastlineWorld, projection=p, col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=putp3p"
mapPlot(coastlineWorld, projection=p, col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=putp5p"
mapPlot(coastlineWorld, projection=p, col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=putp6p"
mapPlot(coastlineWorld, projection=p, col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=qua_aut"
mapPlot(coastlineWorld, projection=p, col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=qsc +lon_0=-100"
mapPlot(coastlineWorld, projection=p, grid=15, col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=robin"
mapPlot(coastlineWorld, projection=p, col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=rouss +lon_0=-40"
mapPlot(coastlineWorld, projection=p, longitudelim=c(-80,0), latitudelim=c(0,45), col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=sinu"
mapPlot(coastlineWorld, projection=p, col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=stere +lat_0=90"
mapPlot(coastlineWorld, projection=p, longitudelim=c(-80,0), latitudelim=c(70, 110), col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=sterea +lat_0=90"
mapPlot(coastlineWorld, projection=p, longitudelim=c(-80,0), latitudelim=c(70, 110), col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=tcea +lon_0=-40"
mapPlot(coastlineWorld, projection=p, longitudelim=c(-80,0), latitudelim=c(0,60), col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=tissot +lat_1=20 +lat_2=60 +lon_0=-40"
mapPlot(coastlineWorld, projection=p, longitudelim=c(-80,0), latitudelim=c(0,60), col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=tmerc +lon_0=-40 +lat_1=30 +lat_2=60"
mapPlot(coastlineWorld, projection=p, longitudelim=c(-80,0), latitudelim=c(0,60), col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=tpeqd +lat_1=30 +lon_1=-80"
mapPlot(coastlineWorld, projection=p, longitudelim=c(-80,0), latitudelim=c(0,60), col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
4
p <- "+proj=tpers +h=1e8"
mapPlot(coastlineWorld, projection=p, col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)
mtext("Problem in Antarctica", line=line, adj=0, col=ecol, font=font)

center

1
2
3
p <- "+proj=ups +ellps=WGS84"
mapPlot(coastlineWorld, projection=p, longitudelim=c(-80,0), latitudelim=c(70, 110), col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=urmfps +n=0.9"
mapPlot(coastlineWorld, projection=p, col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=utm +ellps=WGS84 +lon_0=-40"
mapPlot(coastlineWorld, projection=p, longitudelim=c(-80,0), latitudelim=c(0,45), col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=vandg"
mapPlot(coastlineWorld, projection=p, col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=vitk1 +lat_1=20 +lat_2=60 +lon_0=-40"
mapPlot(coastlineWorld, projection=p, longitudelim=c(-80,0), latitudelim=c(0,45), col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=wag1"
mapPlot(coastlineWorld, projection=p, col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=wag2"
mapPlot(coastlineWorld, projection=p, col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=wag3"
mapPlot(coastlineWorld, projection=p, col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=wag4"
mapPlot(coastlineWorld, projection=p, col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=wag5"
mapPlot(coastlineWorld, projection=p, col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=wag6"
mapPlot(coastlineWorld, projection=p, col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=weren"
mapPlot(coastlineWorld, projection=p, col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=wink1"
mapPlot(coastlineWorld, projection=p, col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=wintri"
mapPlot(coastlineWorld, projection=p, col=col)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

References and resources

  1. Oce website

  2. Jekyll source code for this blog entry: 2015-04-03-oce-proj.Rmd

This website is written in Jekyll, and the source is available on GitHub.