(lat, lon) WKT coordinates do not reproject well with st_transform

Question

I am having some issues importing a file with wkt multipoint features with SRID 4326, for which the coordinates are in order (lat, lon):

>st_crs(4326) 
Coordinate Reference System:
  User input: EPSG:4326 
  wkt:
GEOGCRS["WGS 84",
DATUM["World Geodetic System 1984",
    ELLIPSOID["WGS 84",6378137,298.257223563,
        LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
    ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
    AXIS["geodetic latitude (Lat)",north,
        ORDER[1],
        ANGLEUNIT["degree",0.0174532925199433]],
    AXIS["geodetic longitude (Lon)",east,
        ORDER[2],
        ANGLEUNIT["degree",0.0174532925199433]],
USAGE[
    SCOPE["Horizontal component of 3D system."],
    AREA["World."],
    BBOX[-90,-180,90,180]],
ID["EPSG",4326]]

So I "load" and assign the crs as follows (just one row shown for reproducible example, but this will be hundreds of lines

tst <- data.frame(ID = rep("Test", 2), 
       SRIDTrail = rep(4326, 2), 
       Trail = c("MULTIPOINT (52.86 -8.00, 52.85 -7.98, 52.85 -7.98, 52.85 -7.98, 52.89 -7.98, 52.85 -7.98, 52.85 -7.98, 52.85 -7.98, 52.85 -7.98, 52.85 -7.98, 52.85 -7.98, 52.85 -7.97)", 
                 "MULTIPOINT (53.86 -7.00, 52.02 -6.98, 53.85 -7.80, 51.85 -8.98, 52.89 -7.98, 52.85 -7.98, 52.85 -7.98, 52.85 -7.98, 52.85 -7.98, 52.85 -7.98, 52.85 -7.98, 52.85 -7.97)"))

tst_sf <- tst %>% 
  st_as_sf(wkt = "Trail") %>% 
  st_set_crs(4326)

Now, let's download the world map from the naturalearth package, and check its CRS:

library(rnaturalearth)
world <- ne_countries(scale = "medium", returnclass = "sf")
st_crs(world)

which gives

Coordinate Reference System:
  User input: +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
  wkt:
BOUNDCRS[
    SOURCECRS[
        GEOGCRS["unknown",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ID["EPSG",6326]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8901]],
        CS[ellipsoidal,2],
            AXIS["longitude",east,
                ORDER[1],
                ANGLEUNIT["degree",0.0174532925199433,
                    ID["EPSG",9122]]],
            AXIS["latitude",north,
                ORDER[2],
                ANGLEUNIT["degree",0.0174532925199433,
                    ID["EPSG",9122]]]]],
TARGETCRS[
    GEOGCRS["WGS 84",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        CS[ellipsoidal,2],
            AXIS["latitude",north,
                ORDER[1],
                ANGLEUNIT["degree",0.0174532925199433]],
            AXIS["longitude",east,
                ORDER[2],
                ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4326]]],
ABRIDGEDTRANSFORMATION["Transformation from unknown to WGS84",
    METHOD["Geocentric translations (geog2D domain)",
        ID["EPSG",9603]],
    PARAMETER["X-axis translation",0,
        ID["EPSG",8605]],
    PARAMETER["Y-axis translation",0,
        ID["EPSG",8606]],
    PARAMETER["Z-axis translation",0,
        ID["EPSG",8607]]]]

which states that the first axis corresponds to longitude, not latitude. So, first thing I try (because it'll be faster) is to transform my data to the same projection of the world map, and then plot them both:

tst_sf2 <- tst_sf %>% st_transform(st_crs(world))
ggplot(tst_sf2) + 
  geom_sf(data = world) +
  geom_sf(col = "red") + 
  theme_bw()

This has not worked, as the point that was supposed to be in Ireland is plotted in the Indian ocean, the location with the "swapped" coordinates, that is, lat -8, lon 53).

Let's try the other way round, transform the world map, and not the wtk.

world2 <- world %>% st_transform(st_crs(tst_sf))
ggplot(tst_sf) + 
  geom_sf(data = world2) +
  geom_sf(col = "red") + 
  theme_bw()

This still does not work:

So, my questions are:

(1) Is there any EPSG code I can use that could make R understand the coordinates in the WKT file are swapped with respect to what is expected (I don't mean to enter into the discussion of which order should it be, just to fix it!)

(2) In case that's not possible, how can I change the order of the coordinates, taking into account there will be hundreds of rows and not all multipoints features will be the same length.

Load points as sf, extract coordinates as regular dataframe fields, st_drop_geometry, then load again as sf object with st_as_sf() but with coords in the correct order? E.g., `tst_sf %>% mutate(lat = unlist(map(tst_sf$geometry,1)), lon = unlist(map(tst_sf$geometry,2))) %>% st_drop_geometry() %>% st_as_sf(coords = c("lon", "lat"))` — Skaqqs, Oct 29 '21 at 16:32
@Skaqqs that does not work for multipoints sadly! This only turns it into a single point (POINT feature) where the first coordinate is the latitude of the second point in the multipoint, and the second coordinate is the latitude of the first point in the multipoint, `POINT (52.85631 52.86765)` — Virginia Morera Pujol, Oct 29 '21 at 16:50
for me, the first step, `tst` to `tst_sf` doesn't work: `ogr: corrupt data` — D.J, Oct 29 '21 at 17:04
Convert multipoint to point but retain a grouping variable. Use my above suggestion, then cast point to multipoint based on group as the last step? — Skaqqs, Oct 29 '21 at 17:22
@D.J yeah I realised that later, it's because I split the coordinates in three rows for visibility, if you put them all in the same line it works. — Virginia Morera Pujol, Nov 01 '21 at 09:43

lovalery · Answer 1 · 2021-11-09T20:00:00.580

Please find another solution that takes advantage of the argument authority_compliant = st_axis_order(FALSE/TRUE) of the sf_project() function.

Your data

# The map
world <- ne_countries(scale = "medium", returnclass = "sf")

# Your point(s)
tst <- data.frame(ID = "Test",
                  SRIDTrail = 4326,
                  Trail = "MULTIPOINT (52.86 -8.00, 52.85 -7.98, 52.85 -7.98, 52.85 -7.98, 52.89 -7.98, 52.85 -7.98, 52.85 -7.98, 52.85 -7.98, 52.85 -7.98, 52.85 -7.98, 52.85 -7.98, 52.85 -7.97)")

tst_sf <- tst %>% 
  st_as_sf(wkt = "Trail", crs = 4326) # please, note that I lightened your code a little bit here.

Code

# set the new geometry (i.e. lon-lat instead of lat-long)
RightOrder <- sf_project(from = st_crs(tst_sf), 
                         to = st_crs(world), 
                         matrix(unlist(tst_sf$Trail), 
                                nrow = lapply(tst_sf$Trail, length)[[1]]/2, 
                                ncol = 2), 
                         authority_compliant = st_axis_order(TRUE)) %>% # the argument that allows to choose the order of the axes: lat-lon (FALSE) and lon-lat (TRUE) 
  as.data.frame() %>% 
  setNames(., c("lon", "lat")) %>% 
  st_as_sf(coords = c("lon", "lat"), crs = 4326) %>% 
  st_cast("MULTIPOINT") %>% 
  st_union()

# drop the previous geometry and add the new one
tst_sf <- tst_sf %>%
  st_drop_geometry() %>%
  st_sf(.,RightOrder)

# visualize the result
ggplot(tst_sf) + 
  geom_sf(data = world) +
  geom_sf(col = "red") + 
  theme_bw()

Result

EDIT

Update of the above answer to manage a dataset with multiple rows (cf. comments below)

Please find below the following reprex:

Your data

# The map
world <- ne_countries(scale = "medium", returnclass = "sf")

# Your point(s)
tst <- data.frame(ID = rep("Test", 2), 
                   SRIDTrail = rep(4326, 2), 
                   Trail = c("MULTIPOINT (52.86 -8.00, 52.85 -7.98, 52.85 -7.98, 52.85 -7.98, 52.89 -7.98, 52.85 -7.98, 52.85 -7.98, 52.85 -7.98, 52.85 -7.98, 52.85 -7.98, 52.85 -7.98, 52.85 -7.97)", 
                             "MULTIPOINT (53.86 -7.00, 52.02 -6.98, 53.85 -7.80, 51.85 -8.98, 52.89 -7.98, 52.85 -7.98, 52.85 -7.98, 52.85 -7.98, 52.85 -7.98, 52.85 -7.98, 52.85 -7.98, 52.85 -7.97)"))
                  

tst_sf <- tst %>%  
  st_as_sf(wkt = "Trail", crs = 4326) # please, note that I lightened your code a little bit here.

Code

# set the new geometry (i.e. lon-lat instead of lat-long)
for (i in seq(tst_sf$Trail)){
  
  tst_sf$Trail[i] <- sf_project(from = st_crs(tst_sf), 
                                to = st_crs(world), 
                                matrix(unlist(tst_sf$Trail[i]), 
                                       nrow = lapply(tst_sf$Trail[i], length)[[1]]/2, 
                                       ncol = 2), 
                                authority_compliant = st_axis_order(TRUE)) %>% # the argument that allows to choose the order of the axes: lat-lon (FALSE) and lon-lat (TRUE) 
    as.data.frame() %>% 
    setNames(., c("lon", "lat")) %>% 
    st_as_sf(coords = c("lon", "lat"), crs = 4326) %>% 
    st_cast("MULTIPOINT") %>% 
    st_union()
}

# visualize the result
ggplot(tst_sf) + 
  geom_sf(data = world) +
  geom_sf(col = "red") + 
  theme_bw()

Result

^{Created on 2021-11-09 by the reprex package (v2.0.1)}

This worked but I had to make it a loop/function to run rowwise, because it does not work when the dataset has several rows, each one of them a multipoint! Happy to accept it as an answer if you want to add that, otherwise I'll add an answer myself with the loop when I have a second to tidy up the code! — Virginia Morera Pujol, Nov 02 '21 at 16:50
O.K. I saw that you added some extra rows in your example. So, I'll try to modify the code so that it applies to the different rows. Cheers — lovalery, Nov 02 '21 at 17:04
@Virginia Morera Pujol, I thought you had added MULTIPOINT lines in the data you present in your question. I was wrong! Would it be possible for you to make an edit at the end of your question to bring in a dataset with a few extra rows so that I can do some testing. Ideally, I will try to provide you with a vectorized solution rather than a `for` loop; this will be more efficient. — lovalery, Nov 02 '21 at 17:50
Hello @Virginia Morera Pujol. Thank you very much for the update of your dataset. So, please find my edit at the end of the answer above. I finally opted for a `for` loop because the `project()` function is quite complicated to vectorize. That said, I was able to optimize the code by removing the step that consisted to delete the old column containing the `geometry` and to add the new `geometry` column. To do this, I update the `Trail` column of your `tst_sf` object directly in the `for` loop. Cheers. — lovalery, Nov 09 '21 at 20:04

score 1 · Answer 2 · answered Oct 29 '21 at 20:54

Here is a brute force method, not the most efficient but it seems to work.

tst <- data.frame(ID = "Test", 
                  SRIDTrail = 4326, 
                  Trail = "MULTIPOINT (52.86 -8.00, 52.85 -7.98, 52.85 -7.98, 52.85 -7.98, 52.89 -7.98, 52.85 -7.98, 52.85 -7.98, 52.85 -7.98, 52.85 -7.98, 52.85 -7.98, 52.85 -7.98, 52.85 -7.97)")

tst_sf <- tst %>% st_as_sf(wkt = "Trail") %>% st_set_crs(4326)

#get the latitude & longitude
coordinates <- unlist(tst_sf$Trail)
lat <- coordinates [1:(length(coordinates)/2)]
lon <- coordinates [(length(coordinates)/2+1):length(coordinates)]

#rearrange the columns
#convert back into MULTIPOINT and 
tst_sf$Trail <-sfheaders::sfc_multipoint( matrix(c(lon, lat), ncol=2, byrow=FALSE) )
#redefine the CRS
tst_sf <- tst_sf %>% 
   st_as_sf(wkt = "Trail") %>% 
   st_set_crs(4326)

library(rnaturalearth)
world <- ne_countries(country = 'ireland', scale = "medium", returnclass = "sf")
#st_crs(world)
tst_sf2 <- tst_sf %>% st_transform(st_crs(world))
ggplot(tst_sf2) + 
   geom_sf(data = world) +
   geom_sf(col = "red") + 
   theme_bw()

(lat, lon) WKT coordinates do not reproject well with st_transform

2 Answers2