[This article was first published on Dan Kelley Blog/R, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
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" <- 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, =) |
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, =)
|
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, =)
|
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, =)
mtext("Bad coastlines", line=line, adj=0, col='red')
|
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, =) |
1 2 3 | p <- "+proj=bonne +lat_1=45" mapPlot(coastlineWorld, projection=p) mtext(p, line=line, adj=1, col=pcol, =) |
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, =)
mtext("Bad merdians", line=line, adj=0, col='red')
|
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, =) |
1 2 3 | p <- "+proj=cea" mapPlot(coastlineWorld, projection=p) mtext(p, line=line, adj=1, col=pcol, =) |
1 2 3 | p <- "+proj=collg" mapPlot(coastlineWorld, projection=p) mtext(p, line=line, adj=1, col=pcol, =) |
1 2 3 | p <- "+proj=crast" mapPlot(coastlineWorld, projection=p) mtext(p, line=line, adj=1, col=pcol, =) |
1 2 3 | p <- "+proj=eck1" mapPlot(coastlineWorld, projection=p) mtext(p, line=line, adj=1, col=pcol, =) |
1 2 3 | p <- "+proj=eck2" mapPlot(coastlineWorld, projection=p) mtext(p, line=line, adj=1, col=pcol, =) |
1 2 3 | p <- "+proj=eck3" mapPlot(coastlineWorld, projection=p) mtext(p, line=line, adj=1, col=pcol, =) |
1 2 3 | p <- "+proj=eck4" mapPlot(coastlineWorld, projection=p) mtext(p, line=line, adj=1, col=pcol, =) |
1 2 3 | p <- "+proj=eck5" mapPlot(coastlineWorld, projection=p) mtext(p, line=line, adj=1, col=pcol, =) |
1 2 3 | p <- "+proj=eck6" mapPlot(coastlineWorld, projection=p) mtext(p, line=line, adj=1, col=pcol, =) |
1 2 3 | p <- "+proj=eqc" mapPlot(coastlineWorld, projection=p) mtext(p, line=line, adj=1, col=pcol, =) |
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, =) |
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, =)
mtext("Extraneous lines", line=line, adj=0, col=ecol, =)
|
1 2 3 4 | p <- "+proj=fahey"
mapPlot(coastlineWorld, projection=p)
mtext("+proj=fahey", line=line, adj=1, col=pcol, =)
mtext(p, line=line, adj=1, col=pcol, =)
|
1 2 3 | p <- "+proj=fouc" mapPlot(coastlineWorld, projection=p) mtext(p, line=line, adj=1, col=pcol, =) |
1 2 3 | p <- "+proj=fouc_s" mapPlot(coastlineWorld, projection=p) mtext(p, line=line, adj=1, col=pcol, =) |
1 2 3 | p <- "+proj=gall" mapPlot(coastlineWorld, projection=p) mtext(p, line=line, adj=1, col=pcol, =) |
1 2 3 | p <- "+proj=geos +h=1e8" mapPlot(coastlineWorld, projection=p) mtext(p, line=line, adj=1, col=pcol, =) |
1 2 3 | p <- "+proj=gn_sinu +n=6 +m=3" mapPlot(coastlineWorld, projection=p) mtext(p, line=line, adj=1, col=pcol, =) |
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, =) |
1 2 3 | p <- "+proj=goode" mapPlot(coastlineWorld, projection=p) mtext(p, line=line, adj=1, col=pcol, =) |
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, =)
mtext("Overdrawn coastline", line=line, adj=0, col=ecol, =)
|
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, =)
mtext("Overdrawn coastline", line=line, adj=0, col=ecol, =)
|
1 2 3 | p <- "+proj=hatano" mapPlot(coastlineWorld, projection=p) mtext(p, line=line, adj=1, col=pcol, =) |
1 2 3 | p <- "+proj=healpix" mapPlot(coastlineWorld, projection=p) mtext(p, line=line, adj=1, col=pcol, =) |
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, =)
mtext("Unsure on usage", line=line, adj=0, col=ecol, =)
|
1 2 3 | p <- "+proj=igh" mapPlot(coastlineWorld, projection=p) mtext(p, line=line, adj=1, col=pcol, =) |
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, =)
mtext("Extraneous lines", line=line, adj=0, col=ecol, =)
|
1 2 3 | p <- "+proj=kav5" mapPlot(coastlineWorld, projection=p) mtext(p, line=line, adj=1, col=pcol, =) |
1 2 3 | p <- "+proj=kav7" mapPlot(coastlineWorld, projection=p) mtext(p, line=line, adj=1, col=pcol, =) |
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, =)
mtext("Overdrawn coastline", line=line, adj=0, col=ecol, =)
|
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, =)
|
1 2 3 | p <- "+proj=lonlat"
mapPlot(coastlineWorld, projection=p)
mtext("+proj=lonlat", line=line, adj=1, col=pcol, =)
|
1 2 3 | p <- "+proj=latlon"
mapPlot(coastlineWorld, projection=p)
mtext("+proj=lonlat", line=line, adj=1, col=pcol, =)
|
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, =) |
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, =)
mtext("Extraneous lines", line=line, adj=0, col=ecol, =)
|
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, =) |
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, =)
mtext("Overdrawn coastline", line=line, adj=0, col=ecol, =)
|
1 2 3 | p <- "+proj=loxim" mapPlot(coastlineWorld, projection=p) mtext(p, line=line, adj=1, col=pcol, =) |
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, =)
|
1 2 3 | p <- "+proj=mbt_s" mapPlot(coastlineWorld, projection=p) mtext(p, line=line, adj=1, col=pcol, =) |
1 2 3 | p <- "+proj=mbt_fps" mapPlot(coastlineWorld, projection=p) mtext(p, line=line, adj=1, col=pcol, =) |
1 2 3 | p <- "+proj=mbtfpp" mapPlot(coastlineWorld, projection=p) mtext(p, line=line, adj=1, col=pcol, =) |
1 2 3 | p <- "+proj=mbtfpq" mapPlot(coastlineWorld, projection=p) mtext(p, line=line, adj=1, col=pcol, =) |
1 2 3 | p <- "+proj=mbtfps" mapPlot(coastlineWorld, projection=p) mtext(p, line=line, adj=1, col=pcol, =) |
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, =) |
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, =) |
1 2 3 | p <- "+proj=mill" mapPlot(coastlineWorld, projection=p) mtext(p, line=line, adj=1, col=pcol, =) |
1 2 3 | p <- "+proj=moll" mapPlot(coastlineWorld, projection=p) mtext(p, line=line, adj=1, col=pcol, =) |
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, =) |
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, =)
mtext("Extraneous lines", line=line, adj=0, col=ecol, =)
|
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, =) |
1 2 3 | p <- "+proj=natearth" mapPlot(coastlineWorld, projection=p) mtext(p, line=line, adj=1, col=pcol, =) |
1 2 3 | p <- "+proj=nell" mapPlot(coastlineWorld, projection=p) mtext(p, line=line, adj=1, col=pcol, =) |
1 2 3 | p <- "+proj=nell_h" mapPlot(coastlineWorld, projection=p) mtext(p, line=line, adj=1, col=pcol, =) |
1 2 3 | p <- "+proj=nsper +h=90000000" mapPlot(coastlineWorld, projection=p) mtext(p, line=line, adj=1, col=pcol, =) |
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, =)
mtext("Overdrawn coastline", line=line, adj=0, col=ecol, =)
|
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, =)
mtext("Usage?", line=line, adj=0, col=ecol, =)
|
1 2 3 | p <- "+proj=ocea" mapPlot(coastlineWorld, projection=p) mtext(p, line=line, adj=1, col=pcol, =) |
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, =)
mtext("Extraneous lines", line=line, adj=0, col=ecol, =)
|
1 2 3 | p <- "+proj=ortho" mapPlot(coastlineWorld, projection=p) mtext(p, line=line, adj=1, col=pcol, =) |
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, =)
mtext("Extraneous lines", line=line, adj=0, col=ecol, =)
|
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, =) |
1 2 3 | p <- "+proj=putp1" mapPlot(coastlineWorld, projection=p) mtext(p, line=line, adj=1, col=pcol, =) |
1 2 3 | p <- "+proj=putp2" mapPlot(coastlineWorld, projection=p) mtext(p, line=line, adj=1, col=pcol, =) |
1 2 3 | p <- "+proj=putp3" mapPlot(coastlineWorld, projection=p) mtext(p, line=line, adj=1, col=pcol, =) |
1 2 3 | p <- "+proj=putp5" mapPlot(coastlineWorld, projection=p) mtext(p, line=line, adj=1, col=pcol, =) |
1 2 3 | p <- "+proj=putp6" mapPlot(coastlineWorld, projection=p) mtext(p, line=line, adj=1, col=pcol, =) |
1 2 3 | p <- "+proj=putp3p" mapPlot(coastlineWorld, projection=p) mtext(p, line=line, adj=1, col=pcol, =) |
1 2 3 | p <- "+proj=putp5p" mapPlot(coastlineWorld, projection=p) mtext(p, line=line, adj=1, col=pcol, =) |
1 2 3 | p <- "+proj=putp6p" mapPlot(coastlineWorld, projection=p) mtext(p, line=line, adj=1, col=pcol, =) |
1 2 3 | p <- "+proj=qua_aut" mapPlot(coastlineWorld, projection=p) mtext(p, line=line, adj=1, col=pcol, =) |
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, =) |
1 2 3 | p <- "+proj=robin" mapPlot(coastlineWorld, projection=p) mtext(p, line=line, adj=1, col=pcol, =) |
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, =) |
1 2 3 | p <- "+proj=sinu" mapPlot(coastlineWorld, projection=p) mtext(p, line=line, adj=1, col=pcol, =) |
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, =) |
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, =) |
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, =)
mtext("Overdrawn coastline", line=line, adj=0, col=ecol, =)
|
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, =)
mtext("Extraneous lines", line=line, adj=0, col=ecol, =)
|
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, =) |
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, =) |
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, =) |
1 2 3 | p <- "+proj=tpers +h=1e8" mapPlot(coastlineWorld, projection=p, grid=FALSE) mtext(p, line=line, adj=1, col=pcol, =) |
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, =) |
1 2 3 | p <- "+proj=urmfps +n=0.9" mapPlot(coastlineWorld, projection=p, grid=FALSE) mtext(p, line=line, adj=1, col=pcol, =) |
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, =) |
1 2 3 | p <- "+proj=vandg" mapPlot(coastlineWorld, projection=p, grid=FALSE) mtext(p, line=line, adj=1, col=pcol, =) |
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, =) |
1 2 3 | p <- "+proj=wag1" mapPlot(coastlineWorld, projection=p, grid=FALSE) mtext(p, line=line, adj=1, col=pcol, =) |
1 2 3 | p <- "+proj=wag2" mapPlot(coastlineWorld, projection=p, grid=FALSE) mtext(p, line=line, adj=1, col=pcol, =) |
1 2 3 | p <- "+proj=wag3" mapPlot(coastlineWorld, projection=p, grid=FALSE) mtext(p, line=line, adj=1, col=pcol, =) |
1 2 3 | p <- "+proj=wag4" mapPlot(coastlineWorld, projection=p, grid=FALSE) mtext(p, line=line, adj=1, col=pcol, =) |
1 2 3 | p <- "+proj=wag5" mapPlot(coastlineWorld, projection=p, grid=FALSE) mtext(p, line=line, adj=1, col=pcol, =) |
1 2 3 | p <- "+proj=wag6" mapPlot(coastlineWorld, projection=p, grid=FALSE) mtext(p, line=line, adj=1, col=pcol, =) |
1 2 3 | p <- "+proj=weren" mapPlot(coastlineWorld, projection=p, grid=FALSE) mtext(p, line=line, adj=1, col=pcol, =) |
1 2 3 | p <- "+proj=wink1" mapPlot(coastlineWorld, projection=p, grid=FALSE) mtext(p, line=line, adj=1, col=pcol, =) |
References and resources
-
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 about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.