2

I am trying to define a local CRS in R based on the information in the screenshot image below enter image description here

My R code is as follows:

  1. I define the local CRS from the screenshot image
  2. I apply it to my dummy dataset (1 point)
  3. I apply a conversion back to NAD83

However my output doesn't match. Am I missing information when I do step 1 - or simply doing this wrong?

R Code:

library(sf)
library(rgdal)
library(sp)

df <- data.frame(x_local = 8064.8237, y_local = 1127.1598)
projcrs <- "+proj=tmerc +a=20925604.47417 +ellps=GRS80 +datum=NAD83 +units=us-ft +x_0=32282.466 +y_0=-9739.876 +k_0=1.000115 +lat_0=35.00 +lon_0=-117.65"

sf.df <- st_as_sf(x = df,                         
               coords = c("x_local", "y_local"),
               crs = projcrs)


> new.sf <- st_transform(sf.df, CRS("+proj=longlat +datum=NAD83")@projargs)
> new.sf
Simple feature collection with 1 feature and 0 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: -117.977 ymin: 35.09044 xmax: -117.977 ymax: 35.09044
epsg (SRID):    4269
proj4string:    +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs
                   geometry
1 POINT (-117.977 35.09044)

Whereas the coordinates that I'm expecting are N35°01'47.37315" W117°43'51.16782" (or in decimal degrees 35.0294763 -117.7307693), which is about 25km away!

  • Can you try any more known points? Is it a consistent shift in the same direction? – Spacedman Jul 4 at 15:57
  • @Spacedman: Thanks for asking that. It seems to be a constant shift to the left (WNW) of ~23km. – val Jul 4 at 17:35
3

All the parameters need to be in metres even though your units are survey feet.

See section "False Easting/Northing" in https://proj.org/usage/projections.html where it says "Note that these values are always expressed in metres even if the coordinate system is some other units". Most other parameters do this as well.

I did some conversions...

With:

crs2 <- "+proj=tmerc +a=6378137 +ellps=GRS80 +datum=NAD83 +units=us-ft +x_0=9839.715 +y_0=-2968.72 +k_0=1.000115 +lat_0=35.00 +lon_0=-117.65"

Then:

> sf2.df <- st_as_sf(x = df, coords=1:2, crs=crs2)

Your points go to:

> st_transform(sf2.df, "+init=epsg:4326")
Simple feature collection with 1 feature and 0 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: -117.7309 ymin: 35.02983 xmax: -117.7309 ymax: 35.02983
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
                    geometry
1 POINT (-117.7309 35.02983)

Note that specifying +a and +ellps has a redundancy since the ellipse specifies a semi-major axis itself. In this case +a is the same as the GRS80 semi-major axis so that's not a problem.

I'd check with a few other points to see how this projection works, there may be other things like datum transforms I've not included. But the main principle, that nearly all parameter to PROJ.4 strings are metres, should help.

  • Many thanks! That seems to have been it. I've checked two points so far and they are bang on. Do you know where the documentation on having to list parameters in "meters" comes from? I did not find it when reading the proj.org documentation – val Jul 4 at 17:42
  • 1
    Good point - see edit! – Spacedman Jul 4 at 18:22

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.