Search code examples
rrastercoordinate-systemsutmepsg

How to receive EPSG code from PROJ4 string using R?


How to receive the EPSG code from a PROJ4 string using R?

I have a bunch of raster datasets (randomly distributed all over the world, northern and southern hemisphere), with coordinate reference system UTM WGS84. Now I would like to find out the specific EPSG code for every dataset using R. Do you have any ideas?

# Install and load package raster
install.packages('raster')
library(raster)

# Load a specific raster dataset and query the crs
raster_dat <- raster('D:/UnknownUser/GlobalRasterData/raster_dat1')
crs(raster_dat)
# CRS arguments:
#  +proj=utm +zone=13 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0

Maybe it is possible to create the EPSG code in a dynamic way, since it seems like the EPSG code is always a combination of <326> for the northern hemisphere and the specific zone number e.g. <32> (for Germany) which results in a EPSG code <326> + <32> = <32632>.

Unfortunately i figured out that it seems like you start with <327> for the southern hemisphere and add the specific zone number e.g. <32> (for Nigeria) which results in <327> + <32> =<32732>.

Now I am not sure whether this theory is correct (verification needed incl. source please) and if there is maybe an easier way to receive the EPSG code from the PROJ4 string.

Helpful links I have found:

URL:https://spatialreference.org/ref/epsg/?page=9&search=wgs+84+utm

URL:https://epsg.io/32632

URL:https://de.wikipedia.org/wiki/UTM-Koordinatensystem#/media/Datei:Utm-zones.jpg

EDIT:

# Install and load package rgdal
install.packages('rgdal')
library(rgdal)

# following function could be helpful; it creates a dataframe with 3columns; epsg-code, annotation and proj4-string - unfortunately it seems like there is not an epsg code for each utm wgs84 zone...
epsg_list <- rgdal::make_EPSG()

Thanks a lot in advance, ExploreR


Solution

  • The rgdal package also has a showEPSG() function.

    
    ## example data from sf package
    nc <- sf::read_sf(system.file("shape/nc.shp", package="sf"))
    nc_crs <- raster::crs(nc)
    nc_crs
    #> [1] "+proj=longlat +datum=NAD27 +no_defs"
    
    rgdal::showEPSG(nc_crs)
    #> [1] "4267"
    
    ## Double check it with sf::st_crs
    sf::st_crs(nc)
    #> Coordinate Reference System:
    #>   EPSG: 4267 
    #>   proj4string: "+proj=longlat +datum=NAD27 +no_defs"
    
    ## Your proj4 string
    p4_example <- '+proj=utm +zone=13 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0;'
    rgdal::showEPSG(p4_example)
    #> [1] "32613"
    

    Created on 2020-01-09 by the reprex package (v0.3.0)