r - Setting custom projection in PROJ >=5 paradigm

Question

Before PROJ >=5.0, one could set a custom projection centered anywhere in the world by adjusting the +lon_0=and +lat_0= arguments in proj4string specification.

For example, a projection centered on Europe.

proj <- CRS("+proj=laea +lat_0=53 +lon_0=9 +x_0=4321000 +y_0=3210000 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs")

As I understand it, the use of proj4string is being phased out, and this sort of use is discouraged.

So, with the new projection specification methods in R (e.g. WKT2) is there way to set a custom-centered projection? And if not, is anyone aware of a different automated way of selecting a roughly appropriate projection based on a set of X,Y coordinates (e.g. the centroid)?

Answer 1

You could template the WKT2 version i.e.

lon_0 <- 9
lat_0 <- 53
x_0 <- 4321000
y_0 <- 3210000
template <- "+proj=laea +lat_0={lat_0} +lon_0={lon_0} +x_0={x_0} +y_0={y_0} +datum=WGS84 +units=m +no_defs"

EDIT: I didn't mention dropping towgs84 and ellps here, my take is if WGS84 is the datum and is recognized for later transformation then use that, if towgs84 is all zeroes just drop it. The 2020-campaign against projstrings is ostensibly explained in these posts:

https://www.r-spatial.org/r/2020/03/17/wkt.html

https://cran.r-project.org/web/packages/rgdal/vignettes/PROJ6_GDAL3.html

Now convert that to wkt2 with sf::st_crs(glue::glue(template))$wkt (only since sf Nov 2019 or so), and swap out the parameters with similar text replacement. That way if you need a specific local datum or other authority tweaks you can modify the appropriate wkt2 version.

(sprintf does the same as glue, but the glue way is a bit more convenient)

fwiw, my current understanding is that ellps=WGS84 and towgs84 (with zero params) is unnecessary, 'datum=WGS84' is perfectly well understood by PROJ/GDAL.

Possibly a better starting point, is something like this (I don't know how much the PROJCS name matters if at all). I will need this kind of thing for myself I'm only starting to explore it. I copied this from testepsg (GDAL) output of ESRI:102020.

   lon_0 <- 9
    lat_0 <- 53
    x_0 <- 4321000
    y_0 <- 3210000

tmp <- 'PROJCS["Ad_hoc_handy_wkt2_template_for_laea",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563]],
        PRIMEM["Greenwich",0],
        UNIT["Degree",0.0174532925199433]],
    PROJECTION["Lambert_Azimuthal_Equal_Area"],
    PARAMETER["latitude_of_center",{lat_0}],
    PARAMETER["longitude_of_center",{lon_0}],
    PARAMETER["false_easting",{x_0}],
    PARAMETER["false_northing",{y_0}],
    UNIT["metre",1]]'

sf::st_crs(glue::glue(tmp))
#> Coordinate Reference System:
#>   User input: PROJCS["Ad_hoc_handy_wkt2_template_for_laea",
#> GEOGCS["WGS 84",
#>     DATUM["WGS_1984",
#>         SPHEROID["WGS 84",6378137,298.257223563]],
#>     PRIMEM["Greenwich",0],
#>     UNIT["Degree",0.0174532925199433]],
#> PROJECTION["Lambert_Azimuthal_Equal_Area"],
#> PARAMETER["latitude_of_center",53],
#> PARAMETER["longitude_of_center",9],
#> PARAMETER["false_easting",4321000],
#> PARAMETER["false_northing",3210000],
#> UNIT["metre",1]] 
#>   wkt:
#> PROJCRS["Ad_hoc_handy_wkt2_template_for_laea",
#>     BASEGEOGCRS["WGS 84",
#>         DATUM["World Geodetic System 1984",
#>             ELLIPSOID["WGS 84",6378137,298.257223563,
#>                 LENGTHUNIT["metre",1,
#>                     ID["EPSG",9001]]]],
#>         PRIMEM["Greenwich",0,
#>             ANGLEUNIT["Degree",0.0174532925199433]]],
#>     CONVERSION["unnamed",
#>         METHOD["Lambert Azimuthal Equal Area",
#>             ID["EPSG",9820]],
#>         PARAMETER["Latitude of natural origin",53,
#>             ANGLEUNIT["Degree",0.0174532925199433],
#>             ID["EPSG",8801]],
#>         PARAMETER["Longitude of natural origin",9,
#>             ANGLEUNIT["Degree",0.0174532925199433],
#>             ID["EPSG",8802]],
#>         PARAMETER["False easting",4321000,
#>             LENGTHUNIT["metre",1],
#>             ID["EPSG",8806]],
#>         PARAMETER["False northing",3210000,
#>             LENGTHUNIT["metre",1],
#>             ID["EPSG",8807]]],
#>     CS[Cartesian,2],
#>         AXIS["(E)",east,
#>             ORDER[1],
#>             LENGTHUNIT["metre",1]],
#>         AXIS["(N)",north,
#>             ORDER[2],
#>             LENGTHUNIT["metre",1]]]

^{Created on 2020-06-06 by the reprex package (v0.3.0)}

Caveat, obviously you should check in your context for the appropriateness of a given chunk of metadata like this.