map projections in oce

April 2, 2015
By

(This article was first published on Dan Kelley Blog/R, and kindly contributed to R-bloggers)

Introduction

The latest version (4.9.0) of the PROJ.4 library is being incorporated into the
development version of the oce R package. The work is not finalized yet, but I
thought it might be useful to share an early version of the test suite, so
people could get an idea of the upcoming capabilities.

Note that some projections work quite poorly in oce at the moment. The main
problems are spurious overplotting of coastlines (from the opposite side of the
planet), problems filling land areas that cross longitude cut lines, and
problems filling Antarctica in some projections.

There are many projections to choose from, but in serious work probably only a
few are relevant. Development effort will be focussed on those projections,
and this includes documenting them in oce, so there’s little need to say much
more here.

Note that only PROJ.4 projections that have inverses are incorporated 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 (e.g. the “wintri” projection causes
R to segfault when trying to draw a world coastline).

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
library(oce)
data(coastlineWorld)

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
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)

center

1
2
3
4
p <- "+proj=alsk +lon_0=-150"
mapPlot(coastlineWorld, projection=p, longitudelim=c(-180,-120), latitudelim=c(50,80))
mtext(p, line=line, adj=1, col=pcol, font=font)
mtext("Bad coastlines", line=line, adj=0, col='red')

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("Bad merdians", line=line, adj=0, col='red')

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
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
4
p <- "+proj=gs48"
mapPlot(coastlineWorld, projection=p, longitudelim=c(-130,-60), latitudelim=c(30,60))
mtext(p, line=line, adj=1, col=pcol, font=font)
mtext("Overdrawn coastline", line=line, adj=0, col=ecol, font=font)

center

1
2
3
4
p <- "+proj=gs50"
mapPlot(coastlineWorld, projection=p, longitudelim=c(-130,-60), latitudelim=c(30,60))
mtext(p, line=line, adj=1, col=pcol, font=font)
mtext("Overdrawn coastline", line=line, adj=0, col=ecol, 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
3
4
p <- "+proj=rhealpix"
mapPlot(coastlineWorld, projection=p)# +north_square=1 +south_square=2")
mtext(p, line=line, adj=1, col=pcol, font=font)
mtext("Unsure on usage", 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
4
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)
mtext("Extraneous lines", line=line, adj=0, col=ecol, 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
4
p <- "+proj=krovak +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("Overdrawn coastline", line=line, adj=0, col=ecol, 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
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
4
p <- "+proj=lee_os"
mapPlot(coastlineWorld, projection=p, longitudelim=c(-180,-120), latitudelim=c(70,110))
mtext(p, line=line, adj=1, col=pcol, font=font)
mtext("Overdrawn coastline", line=line, adj=0, col=ecol, 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
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)

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
p <- "+proj=nzmg"
mapPlot(coastlineWorld, projection=p, longitudelim=c(-180,-120), latitudelim=c(-50,-20))
mtext(p, line=line, adj=1, col=pcol, font=font)
mtext("Overdrawn coastline", 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))
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=gstmerc +lon_0=-40 +lat_0=40 +k_0=1"
mapPlot(coastlineWorld, projection=p, longitudelim=c(-80,0), latitudelim=c(0,60))
mtext(p, line=line, adj=1, col=pcol, font=font)
mtext("Overdrawn coastline", line=line, adj=0, col=ecol, 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
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)

center

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

center

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

center

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

center

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

center

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

center

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

center

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

center

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

To leave a comment for the author, please follow the link and comment on their blog: Dan Kelley Blog/R.

R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...



If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Comments are closed.

Search R-bloggers


Sponsors

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)