0

I have selected a custom CRS using the website https://projectionwizard.org/. I chose this custom CRS because I needed a map of the Neotropics (+ or - from the south of the US to the north of Argentina) that preserved distance. Now I have this WKT string and want to add it as the CRS of a raster object inside R. My problem is rather simple, I think: How can I provide R with this WKT in a way it can understand it and create a CRS object?

My custom WKT string is

PROJCS["ProjWiz_Custom_Azimuthal_Equidistant", GEOGCS["GCS_WGS_1984", DATUM["D_WGS_1984", SPHEROID["WGS_1984",6378137.0,298.257223563]], PRIMEM["Greenwich",0.0], UNIT["Degree",0.0174532925199433]], PROJECTION["Azimuthal_Equidistant"], PARAMETER["False_Easting",0.0], PARAMETER["False_Northing",0.0], PARAMETER["Central_Meridian",-82.27], PARAMETER["Latitude_Of_Origin",-12.95], UNIT["Meter",1.0]]

When I do:

sf::st_crs(4326)$wkt

It returns me this:

[1] "GEOGCRS[\"WGS 84\",\n ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n MEMBER[\"World Geodetic System 1984 (Transit)\"],\n MEMBER[\"World Geodetic System 1984 (G730)\"],\n MEMBER[\"World Geodetic System 1984 (G873)\"],\n MEMBER[\"World Geodetic System 1984 (G1150)\"],\n MEMBER[\"World Geodetic System 1984 (G1674)\"],\n MEMBER[\"World Geodetic System 1984 (G1762)\"],\n MEMBER[\"World Geodetic System 1984 (G2139)\"],\n ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1]],\n ENSEMBLEACCURACY[2.0]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n CS[ellipsoidal,2],\n AXIS[\"geodetic latitude (Lat)\",north,\n ORDER[1],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n AXIS[\"geodetic longitude (Lon)\",east,\n ORDER[2],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n USAGE[\n SCOPE[\"Horizontal component of 3D system.\"],\n AREA[\"World.\"],\n BBOX[-90,-180,90,180]],\n ID[\"EPSG\",4326]]"

So I suppose the custom WKT string should be formatted in a way similar to this one (i.e., the WKT string for the WGS84 projection) for R to understand it, but I do not have the faintest idea how can I easily convert it from the format I have. I tried doing it by hand, but I do not understand the syntax very well to do so and it is quite cumbersome anyway.

1 Answer 1

0

You don't need to convert it, its already in a format that st_crs understands.

With your string in an object named proj:

> proj
[1] "PROJCS[\"ProjWiz_Custom_Azimuthal_Equidistant\", GEOGCS[\"GCS_WGS_1984\", DATUM[\"D_WGS_1984\", SPHEROID[\"WGS_1984\",6378137.0,298.257223563]], PRIMEM[\"Greenwich\",0.0], UNIT[\"Degree\",0.0174532925199433]], PROJECTION[\"Azimuthal_Equidistant\"], PARAMETER[\"False_Easting\",0.0], PARAMETER[\"False_Northing\",0.0], PARAMETER[\"Central_Meridian\",-82.27], PARAMETER[\"Latitude_Of_Origin\",-12.95], UNIT[\"Meter\",1.0]]"

I can create a point at some lat-long roughly in the area you are interested in:

> pt = st_point(c(-90, 0))
> pts = st_sfc(pt, crs=4326)

and transform:

> st_transform(pts, st_crs(proj))
Geometry set for 1 feature 
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -867785.2 ymin: 1423362 xmax: -867785.2 ymax: 1423362
Projected CRS: ProjWiz_Custom_Azimuthal_Equidistant
POINT (-867785.2 1423362)

You can also use st_crs(proj)$wkt to see how PROJ has re-interpreted your WKT into a form it seems to prefer.

1

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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