-1

I have two sf objects coming from different sources but supposed to represent the same polygon. However, when plotting them, they clearly appear shifted from one another:

enter image description here

Reprojecting the objects does not seem to solve the problem, and I am wondering whether this is due to one of the object having an "unknown datum based on GRS 1980".

st_crs(sf1)
Coordinate Reference System:
  User input: NAD83 
  wkt:
GEOGCRS["NAD83",
    DATUM["North American Datum 1983",
        ELLIPSOID["GRS 1980",6378137,298.257222101,
            LENGTHUNIT["metre",1]]], [...]


st_crs(sf2)
Coordinate Reference System:
  User input: Unknown datum based upon the GRS 1980 ellipsoid 
  wkt:
GEOGCRS["Unknown datum based upon the GRS 1980 ellipsoid",
    DATUM["Not specified (based on GRS 1980 ellipsoid)",
        ELLIPSOID["GRS 1980",6378137,298.257222101,
            LENGTHUNIT["metre",1]]], [...]

If so, is there any way to identify the datum and make sure that the polygons properly line up?

Reprex with sf objects:

download.file(destfile = "sf1.rds","https://github.com/JoakimWeill/projection_issue_reprex/raw/master/sf1.rds")
download.file(destfile = "sf2.rds","https://github.com/JoakimWeill/projection_issue_reprex/raw/master/sf2.rds")

sf1 <- readRDS("sf1.rds")
sf2 <- readRDS("sf2.rds")

ggplot() +
  geom_sf(data = sf1, fill = "red", alpha = .5) +
  geom_sf(data = sf2, fill = "blue", alpha = .5)

sf1 <- st_transform(sf1, st_crs(sf2))

ggplot() +
  geom_sf(data = sf1, fill = "red", alpha = .5) +
  geom_sf(data = sf2, fill = "blue", alpha = .5)
Matifou
  • 5,399
  • 1
  • 32
  • 41
Jo_
  • 59
  • 5
  • one dirty workaround would have been to check if the difference in X coordinates between sf1 and sf2 is constant, but it looks like the two polygons have a different number of points? `nrow(st_coordinates(sf1))` id 691, and `nrow(st_coordinates(sf2))` is 588? – Matifou Sep 10 '20 at 19:10
  • Yeah these polygons are slightly different (there was some cropping done), I need to find a non-dirty way to figure out if the difference is really due to the different CRS – Jo_ Sep 10 '20 at 19:25

1 Answers1

2

Perhaps you need to define the projection instead of transforming the projection. For example, try defining the datum for the unknown as NAD83 and see if they line up. If not, try defining it as NAD27, transforming to NAD83, then see if it works. You need to tell it the correct starting point before you can transform it into something else. My guess is NAD83 or NAD27 and not WGS84, because of the ellipsoid, and because of how it is shifted. https://www.esri.com/arcgis-blog/products/arcgis-desktop/mapping/wgs84-vs-nad83/

EnKay
  • 48
  • 5
  • The NAD27 datum worked perfectly, thank you! hopefully this will work for all polygons, if not I'll need to find a way to detect the datum from the data – Jo_ Sep 11 '20 at 01:28