2

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()?

4
  • 2
    R should pick up the projection directly from the .prj file with no intervention from you – Ian Turton Mar 12 '20 at 17:10
  • I see it now an attribute of the object Slot "proj4string". – tic-toc-choc Mar 12 '20 at 17:13
  • The standard shapefile coordinate system file suffix is .prj, not .proj – Vince Mar 12 '20 at 17:56
  • I will edit accordingly – tic-toc-choc Mar 12 '20 at 18:50
2

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 %
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))

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.