Find the CRS of a spatial object from .prj

Question

I have two spatial objects (coordinates in .shp files) that I want to draw together in R, it is using the same CRS (i.e. to map these using an identical, appropriate, projection)

(The R code to define a projection is similar to: SpatialPoints(coords, proj4string=CRS('+init=epsg:4326')))

First shape file has an associated .prj that contains: GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]

Second shape file has an associated .prj that contains: PROJCS["WGS_1984_Web_Mercator_Auxiliary_Sphere",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Mercator_Auxiliary_Sphere"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",0.0],PARAMETER["Standard_Parallel_1",0.0],PARAMETER["Auxiliary_Sphere_Type",0.0],UNIT["Meter",1.0]]

I tried to do a search using key terms in this .prj such as https://spatialreference.org/ref/?search=WGS_1984 but don't know which of the several projections to pick.

How can I find the proper argument to CRS()?

Answer 1

They are WKT1:ESRI representations for EPSG:4326 and EPSG:3857 respectively:

C:\>projinfo -o WKT1:ESRI EPSG:4326
WKT1:ESRI string:
GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]


C:\>projinfo -o WKT1:ESRI EPSG:3857
WKT1:ESRI string:
PROJCS["WGS_1984_Web_Mercator_Auxiliary_Sphere",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Mercator_Auxiliary_Sphere"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",0.0],PARAMETER["Standard_Parallel_1",0.0],PARAMETER["Auxiliary_Sphere_Type",0.0],UNIT["Meter",1.0]]

They are very common CRSes, but if you need to find definitions for other WKT representations, you can save them in a text file and get the projinfo identification:

projinfo --identify @the_text_file.txt

Returns:

PROJ.4 string:
+proj=longlat +datum=WGS84 +no_defs +type=crs

WKT2:2019 string:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]],
        ID["EPSG",6326]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["Degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["longitude",east,
            ORDER[1],
            ANGLEUNIT["Degree",0.0174532925199433]],
        AXIS["latitude",north,
            ORDER[2],
            ANGLEUNIT["Degree",0.0174532925199433]]]

Identification match count: 1
EPSG:4326: 100 %

and:

PROJ.4 string:
+proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs +type=crs

WKT2:2019 string:
PROJCRS["WGS 84 / Pseudo-Mercator",
    BASEGEOGCRS["WGS 84",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ID["EPSG",6326]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["Degree",0.0174532925199433]]],
    CONVERSION["unnamed",
        METHOD["Popular Visualisation Pseudo Mercator",
            ID["EPSG",1024]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",0,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["False easting",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]

Identification match count: 1
EPSG:3857: 100 %

Answer 2

Realised thanks to a comment that the projection is already there. I picked one of the two (guyana_regions@proj4string@projargs) and proceed to convert the other smoothly.

library(rgdal)

guyana_regions_shp <- './Data/Guyana/Shapefiles/Region.shp'
guyana_regions <- readOGR(guyana_regions_shp)

suriname_resorts_shp <- './Data/Suriname/Shapefiles/Resorts/Suriname_Ressorts.shp'
suriname_resorts <- readOGR(suriname_resorts_shp)

suriname_resorts <- spTransform(suriname_resorts, CRS(guyana_regions@proj4string@projargs))