map projections in rgdal and sf

Posted on August 1, 2020 by Dan Kelley Blog/R in R bloggers | 0 Comments

[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.

Executive Summary. a transition of oce to sf will require the removal of healpix, nesper, pconic and rhealpix projections, all of which cause errors in sf. The loss of healpix, rhealpix and pconic projections is not of much concern, since these were interrupted projections that were not handled properly by oce::mapPlot() anyway (see http://dankelley.github.io/r/2020/01/04/oce-proj.html for the healpix and rhealpix examples; the nesper case was ignored in that blog item). The nesper projection seems unlikely to enjoy wide usage, since it is showing a near-earth satellite view. Thus, the loss of these 4 oce projections should not stand in the way of a transition from rgdal to sf \ldots and the benefits of switching to the newer sf package are convincing, in terms of accuracy and future availability.

Introduction

For several years, oce has used the rgdal::transform() function for map-projection calculations. This cannot continue, given that the rgdal developers intend to deprecate transform(). Lately, I have been exploring the possibility of switching to sf::st_transform() for map-projection calculations. This work is being done in the sf branch of oce. In the present initial phase of the transition, calculations are done with both rgdal::transform() and sf::st_transform(), and warnings are issued if the results differ. (This requires version 0.8.1 of the st package, which is presently only available on github.)

This blog posting tests the sf branch of oce against the projections used in the 2015-04-03-oce-proj posting. My quick scan shows that this works as expected.

Note that this version of oce accepts lonlat and longlat as synonyms, and the same for latlon and latlong.

library(oce)
data(coastlineWorld)
cl45 <- coastlineCut(coastlineWorld, lon_0=-45)

par(mar=rep(2, 4))
line <- 0.25
pcol <- "blue"
ecol <- "red"
font <- 1
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

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

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

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

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

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

center

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

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

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

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

center

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

center

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

center

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

center

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

center

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

center

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

center

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

center

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

center

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

center

p <- "+proj=eqearth"
mapPlot(coastlineWorld, projection=p, col=col)

## Warning in CPL_crs_parameters(x): GDAL Error 1: PROJ: proj_as_wkt: Unsupported
## conversion method: Equal Earth
## Warning in CPL_crs_parameters(x): GDAL Error 1: PROJ: proj_as_wkt: Unsupported
## conversion method: Equal Earth
mtext(p, line=line, adj=1, col=pcol, font=font)

center

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

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

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

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

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

center

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

center

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

center

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

center

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

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

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

center

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

center

## fails 2020-08-02 ## p <- "+proj=healpix +a=1"
## fails 2020-08-02 ## mapPlot(coastlineWorld, projection=p, col=col)
## fails 2020-08-02 ## mtext(p, line=line, adj=1, col=pcol, font=font)
## fails 2020-08-02 ## mtext("Problem in Canada", line=line, adj=0, col=ecol, font=font)

## fails 2020-08-02 ## p <- "+proj=rhealpix +south_square=0 +north_square=1"
## fails 2020-08-02 ## mapPlot(coastlineWorld, projection=p)
## fails 2020-08-02 ## mtext(p, line=line, adj=1, col=pcol, font=font)
## fails 2020-08-02 ## mtext("Is this an acid trip?", line=line, adj=0, col=ecol, font=font)

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

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

center

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

center

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

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

center

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

center

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

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

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

center

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

center

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

center

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

center

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

center

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

center

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

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

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

center

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

center

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

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

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

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

center

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

center

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

center

## fails 2020-08-02 ## p <- "+proj=nsper +h=90000000"
## fails 2020-08-02 ## mapPlot(coastlineWorld, projection=p, col=col)
## fails 2020-08-02 ## mtext(p, line=line, adj=1, col=pcol, font=font)
## fails 2020-08-02 ## mtext("Problem in Antarctica", line=line, adj=0, col=ecol, font=font)

p <- "+proj=ob_tran +o_proj=mill +o_lon_p=40 +o_lat_p=50"
mapPlot(coastlineWorld, projection=p, col=col)
mtext("Ugly Horizontal Lines; weird gradicules", line=line, adj=0, col=ecol, font=font)

center

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

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

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

## fails 2020-08-02 ## > p <- "+proj=pconic +lat_1=20 +lat_2=60 +lon_0=-40"
## fails 2020-08-02 ## > mapPlot(coastlineWorld, projection=p, longitudelim=c(-80,0), latitudelim=c(0,60), col=col)
## fails 2020-08-02 ## Error in CPL_geos_op2(op, x, y) : 
## fails 2020-08-02 ##   Evaluation error: TopologyException: found non-noded intersection between LINESTRING (2.90185e+06 6.38497e+07, 3.59738e+06 7.62361e+07) and LINESTRING (2.42138e+07 4.43386e+08, -699901 -292340) at 2901848.262176008 63849691.953739122.

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

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

center

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

center

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

center

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

center

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

center

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

center

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

center

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

center

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

center

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

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

center

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

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

center

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

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

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

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

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

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

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

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

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

center

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

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

center

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

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

center

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

center

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

center

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

center

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

center

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

center

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

center

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

center

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. Previous related posting dated 2020-01-04

  3. Jekyll source code for this blog entry: 2020-01-04-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.

Copyright © 2022 | MH Corporate basic by MH Themes

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)