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 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, available on the CRAN archive. However, a CRAN member suggested removing the PROJ.4 code from oce when 0.9-17 was being vetted, and so I started exploring the use of the rgdal package. At the moment, this work is only on the git branch named project_with_rgdal. The results with this version of oce seem to be very close to those with the previous system, with differences only in the cases where extraneous lines are seen.

Note that only PROJ.4 projections that have inverses are available in oce. 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 serives 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
library(oce)
data(coastlineWorld)
cl45 <- coastlineCut(coastlineWorld, lon_0=-45)

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

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

center

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

center

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

center

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

center

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

center

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

center

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

center

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

center

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

center

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

center

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

center

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

center

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

center

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

center

1
2
3
p <- "+proj=eqc"
mapPlot(coastlineWorld, projection=p)
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))
mtext(p, line=line, adj=1, col=pcol, font=font)

center

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

center

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

center

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

center

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

center

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

center

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

center

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

center

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

center

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

center

1
2
p <- "+proj=rhealpix"
mapPlot(coastlineWorld, projection=p)# +north_square=1 +south_square=2")
## Error in seq.default(-90 + small, 90 - small, dlatitude): invalid (to - from)/by in seq(.)
1
2
mtext(p, line=line, adj=1, col=pcol, font=font)
mtext("Extraneous lines", line=line, adj=0, col=ecol, font=font)

center

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

center

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

center

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

center

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

center

1
2
3
p <- "+proj=laea +lon_0=-40"
mapPlot(coastlineWorld, projection=p)
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)
mtext("+proj=lonlat", line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=latlon"
mapPlot(coastlineWorld, projection=p)
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))
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
4
p <- "+proj=lcca +lat_0=20 +lon_0=-40"
mapPlot(coastlineWorld, projection=p, longitudelim=c(-80, 0), latitudelim=c(0, 60))
mtext(p, line=line, adj=1, col=pcol, font=font)
mtext("Extraneous lines", line=line, adj=0, col=ecol, font=font)

center

1
2
3
4
5
p <- "+proj=lcca +lat_0=20 +lon_0=-40"
mapPlot(cl45, projection=p, longitudelim=c(-80, 0), latitudelim=c(0, 60))
mtext(p, line=line, adj=1, col=pcol, font=font)
mtext("Using coastlineCut()", line=line+0.9, adj=1, col=pcol, font=font)
mtext("Extraneous lines", line=line, adj=0, col=ecol, font=font)

center

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

center

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

center

1
2
3
4
5
p <- "+proj=lsat +lsat=1 +path=200"
mapPlot(coastlineWorld, projection=p,
            longitudelim=c(-180,-120), latitudelim=c(70,120))
mtext(p, line=line, adj=1, col=pcol, font=font)
mtext("Extraneous lines", line=line, adj=0, col=ecol, font=font)

center

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

center

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

center

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

center

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

center

1
2
3
p <- "+proj=mbtfps"
mapPlot(coastlineWorld, projection=p)
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))
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))
mtext(p, line=line, adj=1, col=pcol, font=font)

center

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

center

1
2
3
p <- "+proj=moll"
mapPlot(coastlineWorld, projection=p)
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))
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
4
p <- "+proj=murd2 +lat_1=30 +lat_2=60 +lon_0=-40"
mapPlot(coastlineWorld, projection=p, longitudelim=c(-80,0), latitudelim=c(0,45))
mtext(p, line=line, adj=1, col=pcol, font=font)
mtext("Extraneous lines", line=line, adj=0, col=ecol, 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))
mtext(p, line=line, adj=1, col=pcol, font=font)

center

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

center

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

center

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

center

1
2
3
p <- "+proj=nsper +h=90000000"
mapPlot(coastlineWorld, projection=p)
mtext(p, line=line, adj=1, col=pcol, 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))
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
p <- "+proj=ocea"
mapPlot(coastlineWorld, projection=p)
mtext(p, line=line, adj=1, col=pcol, font=font)

center

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

center

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

center

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

center

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

center

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

center

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

center

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

center

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

center

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

center

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

center

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

center

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

center

1
2
3
p <- "+proj=qua_aut"
mapPlot(coastlineWorld, projection=p)
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, fill='lightgray')
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=robin"
mapPlot(coastlineWorld, projection=p)
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))
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
p <- "+proj=sinu"
mapPlot(coastlineWorld, projection=p)
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))
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))
mtext(p, line=line, adj=1, col=pcol, font=font)

center

1
2
3
4
p <- "+proj=tcea +lon_0=-40"
mapPlot(coastlineWorld, projection=p, longitudelim=c(-80,0), latitudelim=c(0,60))
mtext(p, line=line, adj=1, col=pcol, font=font)
mtext("Extraneous lines", line=line, adj=0, col=ecol, 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))
mtext(p, line=line, adj=1, col=pcol, font=font)

center

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

center

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

center

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

center

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

center

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

center

1
2
3
p <- "+proj=vandg"
mapPlot(coastlineWorld, projection=p)
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))
mtext(p, line=line, adj=1, col=pcol, font=font)

center

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

center

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

center

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

center

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

center

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

center

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

center

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

center

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

center

1
2
3
p <- "+proj=wintri"
mapPlot(coastlineWorld, projection=p)
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.