Sf: SRID = 0 when writing WGS 84 shp file to Postgres
When reading in a shp file, which is supposed to have a EPSG 4326 CRS, and writing that to Postgres, the SRID in Postgres is not set correctly. Unfortunately, I cannot share the shp file in question, so I will have to try to reproduce this with the nc.shp.
First, what I would like to get is this:
library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
library(RPostgres)
con <-
dbConnect(
drv = RPostgres::Postgres(),
dbname = "my_db",
host = "my_host",
password = "my_pass",
user = "my_user"
)
nc = st_read(system.file("shape/nc.shp", package="sf"))
#> Reading layer `nc' from data source `C:\Users\daniel\Documents\.R\win-library\sf\shape\nc.shp' using driver `ESRI Shapefile'
#> Simple feature collection with 100 features and 14 fields
#> geometry type: MULTIPOLYGON
#> dimension: XY
#> bbox: xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
#> geographic CRS: NAD27
nc_4326 = st_transform(nc, crs = 4326)
st_crs(nc_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["unknown"],
#> AREA["World"],
#> BBOX[-90,-180,90,180]],
#> ID["EPSG",4326]]
dbWriteTable(
conn = con,
name = "nc_4326",
value = nc_4326,
overwrite = TRUE,
temporary = TRUE
)
#> Note: method with signature 'DBIObject#sf' chosen for function 'dbDataType',
#> target signature 'PqConnection#sf'.
#> "PqConnection#ANY" would also be valid
dbGetQuery(con, statement = "select st_srid(geometry) from nc_4326 limit 1;")
#> st_srid
#> 1 4326
I have dput() the crs of my shp-file and transform the nc data to the this crs.
my_crs <-
structure(
list(
input = "WGS 84",
wkt = "GEOGCRS[\"WGS 84\",\n DATUM[\"World Geodetic System 1984\",\n ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1]]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n CS[ellipsoidal,2],\n AXIS[\"latitude\",north,\n ORDER[1],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n AXIS[\"longitude\",east,\n ORDER[2],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n ID[\"EPSG\",4326]]"),
class = "crs"
)
nc_wgs84 <- st_transform(nc, crs = my_crs)
st_crs(nc_wgs84)
#> Coordinate Reference System:
#> User input: WGS 84
#> 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["latitude",north,
#> ORDER[1],
#> ANGLEUNIT["degree",0.0174532925199433]],
#> AXIS["longitude",east,
#> ORDER[2],
#> ANGLEUNIT["degree",0.0174532925199433]],
#> ID["EPSG",4326]]
There are only slight, but apparently crucial, differences in the crs AFAICT, i.e.聽the User input, AXIS, and USAGE statements differ.
But nothing changes, if I try to transform this to EPSG:4326, so apparently sf (PROJ?) considers them to be equivalent.
st_crs(st_transform(nc_wgs84, crs = 4326))
#> Coordinate Reference System:
#> User input: WGS 84
#> 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["latitude",north,
#> ORDER[1],
#> ANGLEUNIT["degree",0.0174532925199433]],
#> AXIS["longitude",east,
#> ORDER[2],
#> ANGLEUNIT["degree",0.0174532925199433]],
#> ID["EPSG",4326]]
But if I write them to Postgres, the SRID is not written correctly AFAICT.
dbWriteTable(
conn = con,
name = "nc_wgs84",
value = nc_wgs84,
overwrite = TRUE,
temporary = TRUE
)
dbGetQuery(con, statement = "select st_srid(geometry) from nc_wgs84 limit 1;")
#> st_srid
#> 1 0
dbDisconnect(con)
Is this a bug or am I missing something? And if it is a bug, is it in the import (CRS not read correctly from the shp) or the export (writing to Postgres) step?
Given that st_transform(crs = 4326) does not seem to actually transform, forcing the st_crs(4326) ought to be a safe workaround, I suppose?
st_crs(nc_wgs84) <- st_crs(4326)
FWIW, I also tried with a shp file that used to work with earlier (pre v0.9?) versions of sf and I am having the same problem with them now.
Session info
devtools::session_info()
#> - Session info ---------------------------------------------------------------
#> setting value
#> version R version 4.0.2 (2020-06-22)
#> os Windows 10 x64
#> system x86_64, mingw32
#> ui RTerm
#> language en
#> collate German_Germany.1252
#> ctype German_Germany.1252
#> tz Europe/Berlin
#> date 2020-09-17
#>
#> - Packages -------------------------------------------------------------------
#> package * version date lib source
#> assertthat 0.2.1 2019-03-21 [1] CRAN (R 4.0.2)
#> backports 1.1.9 2020-08-24 [1] CRAN (R 4.0.2)
#> bit 4.0.4 2020-08-04 [1] CRAN (R 4.0.2)
#> bit64 4.0.5 2020-08-30 [1] CRAN (R 4.0.2)
#> blob 1.2.1 2020-01-20 [1] CRAN (R 4.0.2)
#> callr 3.4.4 2020-09-07 [1] CRAN (R 4.0.2)
#> class 7.3-17 2020-04-26 [2] CRAN (R 4.0.2)
#> classInt 0.4-3 2020-04-07 [1] CRAN (R 4.0.2)
#> cli 2.0.2 2020-02-28 [1] CRAN (R 4.0.2)
#> crayon 1.3.4 2017-09-16 [1] CRAN (R 4.0.2)
#> DBI 1.1.0 2019-12-15 [1] CRAN (R 4.0.2)
#> desc 1.2.0 2018-05-01 [1] CRAN (R 4.0.2)
#> devtools 2.3.1 2020-07-21 [1] CRAN (R 4.0.2)
#> digest 0.6.25 2020-02-23 [1] CRAN (R 4.0.2)
#> dplyr 1.0.2 2020-08-18 [1] CRAN (R 4.0.2)
#> e1071 1.7-3 2019-11-26 [1] CRAN (R 4.0.2)
#> ellipsis 0.3.1 2020-05-15 [1] CRAN (R 4.0.2)
#> evaluate 0.14 2019-05-28 [1] CRAN (R 4.0.2)
#> fansi 0.4.1 2020-01-08 [1] CRAN (R 4.0.2)
#> fs 1.5.0 2020-07-31 [1] CRAN (R 4.0.2)
#> generics 0.0.2 2018-11-29 [1] CRAN (R 4.0.2)
#> glue 1.4.2 2020-08-27 [1] CRAN (R 4.0.2)
#> highr 0.8 2019-03-20 [1] CRAN (R 4.0.2)
#> hms 0.5.3 2020-01-08 [1] CRAN (R 4.0.2)
#> htmltools 0.5.0 2020-06-16 [1] CRAN (R 4.0.2)
#> KernSmooth 2.23-17 2020-04-26 [2] CRAN (R 4.0.2)
#> knitr 1.29 2020-06-23 [1] CRAN (R 4.0.2)
#> lifecycle 0.2.0 2020-03-06 [1] CRAN (R 4.0.2)
#> magrittr 1.5.0.9000 2020-08-31 [1] Github (tidyverse/magrittr@15f6f07)
#> memoise 1.1.0 2017-04-21 [1] CRAN (R 4.0.2)
#> pillar 1.4.6 2020-07-10 [1] CRAN (R 4.0.2)
#> pkgbuild 1.1.0 2020-07-13 [1] CRAN (R 4.0.2)
#> pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 4.0.2)
#> pkgload 1.1.0 2020-05-29 [1] CRAN (R 4.0.2)
#> prettyunits 1.1.1 2020-01-24 [1] CRAN (R 4.0.2)
#> processx 3.4.4 2020-09-03 [1] CRAN (R 4.0.2)
#> ps 1.3.4 2020-08-11 [1] CRAN (R 4.0.2)
#> purrr 0.3.4 2020-04-17 [1] CRAN (R 4.0.2)
#> R6 2.4.1 2019-11-12 [1] CRAN (R 4.0.2)
#> Rcpp 1.0.5 2020-07-06 [1] CRAN (R 4.0.2)
#> remotes 2.2.0 2020-07-21 [1] CRAN (R 4.0.2)
#> rlang 0.4.7 2020-07-09 [1] CRAN (R 4.0.2)
#> rmarkdown 2.3 2020-06-18 [1] CRAN (R 4.0.2)
#> RPostgres * 1.2.0 2019-12-18 [1] CRAN (R 4.0.2)
#> rprojroot 1.3-2 2018-01-03 [1] CRAN (R 4.0.2)
#> sessioninfo 1.1.1 2018-11-05 [1] CRAN (R 4.0.2)
#> sf * 0.9-6 2020-09-13 [1] CRAN (R 4.0.2)
#> stringi 1.5.3 2020-09-09 [1] CRAN (R 4.0.2)
#> stringr 1.4.0 2019-02-10 [1] CRAN (R 4.0.2)
#> testthat 2.3.2 2020-03-02 [1] CRAN (R 4.0.2)
#> tibble 3.0.3 2020-07-10 [1] CRAN (R 4.0.2)
#> tidyselect 1.1.0 2020-05-11 [1] CRAN (R 4.0.2)
#> units 0.6-7 2020-06-13 [1] CRAN (R 4.0.2)
#> usethis 1.6.1 2020-04-29 [1] CRAN (R 4.0.2)
#> vctrs 0.3.4 2020-08-29 [1] CRAN (R 4.0.2)
#> withr 2.2.0 2020-04-20 [1] CRAN (R 4.0.2)
#> xfun 0.17 2020-09-09 [1] CRAN (R 4.0.2)
#> yaml 2.2.1 2020-02-01 [1] CRAN (R 4.0.0)
#>
#> [1] C:/Users/daniel/Documents/.R/win-library
#> [2] C:/Program Files/R/R-4.0.2/library
dpprdan
All 18 comments
Which postgis version is this?
@etiennebr
edzer
on 17 Sep 2020
Thanks for the report @dpprdan, and thanks for tagging me @edzer. My first thought is that it could be an issue with Windows (which I haven't tested much), or with conflicting versions. As @edzer mentionned, in addition to your session_info(), could you provide:
DBI::dbGetQuery(con, "select postgis_version()")
DBI::dbGetQuery(con, "select version()")
sf::sf_extSoftVersion()
etiennebr
on 17 Sep 2020
I ran your code. You are right, the differences in the wkt create the issue (the synchronization of crs uses the wkt in priority). So it's likely the version of the creator that's different from the one you're running. It would be nice to have a warning, still.
You probably figured a workaround (manually force the crs), but just in case:
dbWriteTable(
conn = con,
name = "nc_wgs84",
value = st_set_crs(nc_wgs84, 4326), # you must be certain that the srid is correct
overwrite = TRUE,
temporary = TRUE
)
Thanks for the swift response! :rocket:
DBI::dbGetQuery(con, "select postgis_version()")
#> postgis_version
#> 1 2.5 USE_GEOS=1 USE_PROJ=1 USE_STATS=1
DBI::dbGetQuery(con, "select version()")
#> version
#> 1 PostgreSQL 11.4 (Ubuntu 11.4-1.pgdg18.10+1) on x86_64-pc-linux-gnu, compiled by gcc (Ubuntu 8.3.0-6ubuntu1~18.10.1) 8.3.0, 64-bit
sf::sf_extSoftVersion()
#> GEOS GDAL proj.4 GDAL_with_GEOS USE_PROJ_H
#> "3.8.0" "3.0.4" "6.3.1" "true" "true"
馃 I thought Postgis would be newer... I saw "We don't test sf against postgis 2.5 anymore". I can try against the docker image tomorrow, probably.
You probably figured a workaround (manually force the crs)
Yes, similar (st_crs(nc_wgs84) <- st_crs(4326)) but yours is more concise. The "be certain" part might be tricky 馃槈. I'd do a st_transform(crs=4326) before, just to be on the safe side.
(the synchronization of crs uses the wkt in priority)
Not sure, I fully understand that, esp. "in priority" (not sure I need to, either).
dpprdan
on 17 Sep 2020
Thanks, I checked against a more recent version, so I'm pretty sure it's not the source of the problem.
Since sf version 0.9 the crs are stored in wkt format. It's the right thing to do because it's more robust to discrepancies with the proj database. However, the interface with postgis hasn't evolved, so we can't just send a srid anymore. We first scan all the database projections and try to match a wkt. If we have a match, then we use it. The database could have a different srid than the local proj database. For that reason, it's better to use the wkt. If we can't find a matching wkt, we'll use the crs IF it was specified by the user (e.g. with st_crs or st_transform) and if the local wkt matches the database wkt. Otherwise, we'll just use 0 (no srid), which is the part that would need a warning. This really should be written somewhere :) I've been planning to work on this since a while, so I won't give a deadline, but I'd like to make it more robust with better crs matching.
etiennebr
on 17 Sep 2020
Thanks @etiennebr! I think I understand it better now.
Apparently st_transform(nc_wgs84, crs = 4326) does not change the wkt string my case 鈽濓笍. I am wondering whether it should? In other words:
If we can't find a matching
wkt, we'll use the crs IF it was specified by the user (e.g. withst_crsorst_transform) and if the localwktmatches the databasewkt.
Both ifs should be matched in the following, should they not?
nc = st_read(system.file("shape/nc.shp", package="sf"))
my_crs <-
structure(
list(
input = "WGS 84",
wkt = "GEOGCRS[\"WGS 84\",\n DATUM[\"World Geodetic System 1984\",\n ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1]]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n CS[ellipsoidal,2],\n AXIS[\"latitude\",north,\n ORDER[1],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n AXIS[\"longitude\",east,\n ORDER[2],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n ID[\"EPSG\",4326]]"),
class = "crs"
)
nc_wgs84 <- st_transform(nc, crs = my_crs)
Until here I have created the sf object with "my_crs". Now I transform to EPSG:4326.
nc_4326 <- st_transform(nc_wgs84, crs = 4326)
But the wkt does not change.
st_crs(nc_4326)
#> Coordinate Reference System:
#> User input: WGS 84
#> 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["latitude",north,
#> ORDER[1],
#> ANGLEUNIT["degree",0.0174532925199433]],
#> AXIS["longitude",east,
#> ORDER[2],
#> ANGLEUNIT["degree",0.0174532925199433]],
#> ID["EPSG",4326]]
And therefore the SRID is not set in Postgis.
dbWriteTable(
conn = con,
name = "nc_4326",
value = nc_4326,
overwrite = TRUE,
temporary = TRUE
)
#> Note: method with signature 'DBIObject#sf' chosen for function 'dbDataType',
#> target signature 'PqConnection#sf'.
#> "PqConnection#ANY" would also be valid
dbGetQuery(con, statement = "select st_srid(geometry) from nc_4326 limit 1;")
#> st_srid
#> 1 0
So the CRS was apparently not used even though I specified it. Possibly because st_transform did not transform, because it considered the CRS to be the same? What does st_transform() match on, i.e. how does it determine whether a transformation is necessary?
Or is this part of "I'd like to make it more robust with better crs matching."?
dpprdan
on 18 Sep 2020
Thanks for your persistence @dpprdan! I'm pretty sure this is a bug.
I'm pretty sure there is something missing in the conversion to binary for the database.
https://github.com/r-spatial/sf/blob/f8e6d27f3deeab543d57b1334b548dc62a6190f2/R/db.R#L500-L502
@edzer, I can't remember the intention behind 13f1c9f where the srid argument was removed. I can look into it, probably next week.
etiennebr
on 18 Sep 2020
Before sf adopted wkt as the item that describes the CRS, we had a crs object that consisted of an srid integer (optional, possibly NA) and proj4string. We abandoned that, and now allow for all kind of input (field input) but a wkt field that is the key. When input is e.g. EPSG:4326, sf will try to recover the number 4326 from it, and use that as srid, which is no longer in the crs; this was done in the section starting here: https://github.com/r-spatial/sf/commit/9f9418318953add34a9741aaab1e67023a2ef29b#diff-7d5aa9bb64a00de8cb92bc8b1e2cdd24R682
edzer
on 18 Sep 2020
I came up with a minimal reprex. The issue is with st_as_binary, which removes any crs that are not in the local database. We used to be able to force this with the srid argument, but it's not there anymore.
#
library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 7.0.0
my_crs <-
structure(
list(
input = "WGS 84",
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["latitude",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["longitude",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4326]]'),
class = "crs"
)
x <- st_as_sfc("POINT(0 0)")
x <- st_set_crs(x, my_crs)
st_crs(x)
#> Coordinate Reference System:
#> User input: WGS 84
#> 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["latitude",north,
#> ORDER[1],
#> ANGLEUNIT["degree",0.0174532925199433]],
#> AXIS["longitude",east,
#> ORDER[2],
#> ANGLEUNIT["degree",0.0174532925199433]],
#> ID["EPSG",4326]]
binx <- st_as_binary(x)
st_crs(st_as_sfc(binx))
#> Coordinate Reference System: NA
Created on 2020-09-19 by the reprex package (v0.3.0)
Probably here...
https://github.com/r-spatial/sf/blob/f8e6d27f3deeab543d57b1334b548dc62a6190f2/src/wkb.cpp#L680-L691
etiennebr
on 19 Sep 2020
@edzer, shouldn't we check the wkt against the local database for a srid before using the input?
proj4 strings also can't be used in st_as_binary().
x <- st_set_crs(x, "+proj=longlat +datum=WGS84 +no_defs")
st_crs(x)
#> Coordinate Reference System:
#> User input: +proj=longlat +datum=WGS84 +no_defs
#> wkt:
#> 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]]]]
binx <- st_as_binary(x)
st_crs(st_as_sfc(binx))
#> Coordinate Reference System: NA
I see two elements:
st_as_binaryshould probably use thewktto find thesrid- The
sridoverride was useful for locally missing projections. I'm surprised none of the tests failed (apparently it wasn't covered).
etiennebr
on 19 Sep 2020
Is it conceivable that part of the trouble is axis order? EPSG:4326 mandates Latitude-Longitude, so will not be seen as such if the axis order is visualization/GIS? Or is it just me, seeing axis-order issues behind any lamp-post?
rsbivand
on 20 Sep 2020
@etiennebr you need to use EWKB=TRUE to roundtrip the CRS:
library(sf)
# Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 7.0.0
x = st_sf(x = 1, geom = st_sfc(st_point(c(0,0)), crs = 4326))
st_as_sfc(st_as_binary(st_geometry(x), EWKB=TRUE), EWKB=TRUE)
# Geometry set for 1 feature
# geometry type: POINT
# dimension: XY
# bbox: xmin: 0 ymin: 0 xmax: 0 ymax: 0
# geographic CRS: WGS 84
# POINT (0 0)
And:
> st_as_sfc(st_as_binary(st_geometry(x), EWKB=TRUE), EWKB=TRUE) %>% st_crs() %>% `$`("epsg")
[1] 4326
Yes, sorry! Here's a better reprex:
library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 7.0.0
x = st_as_sfc("POINT(0 0)", crs = 4326)
st_as_sfc(st_as_binary(x, EWKB=TRUE), EWKB=TRUE)
#> Geometry set for 1 feature
#> geometry type: POINT
#> dimension: XY
#> bbox: xmin: 0 ymin: 0 xmax: 0 ymax: 0
#> geographic CRS: WGS 84
#> POINT (0 0)
x = st_as_sfc("POINT(0 0)", crs = "+proj=longlat +datum=WGS84 +no_defs")
st_as_sfc(st_as_binary(x, EWKB=TRUE), EWKB=TRUE)
#> Geometry set for 1 feature
#> geometry type: POINT
#> dimension: XY
#> bbox: xmin: 0 ymin: 0 xmax: 0 ymax: 0
#> CRS: NA
#> POINT (0 0)
Created on 2020-09-20 by the reprex package (v0.3.0)
etiennebr
on 21 Sep 2020
Yes, short for
> st_crs("+proj=longlat +datum=WGS84 +no_defs")$epsg
[1] NA
or even
> st_crs(st_crs(4326)$proj4string)$epsg
[1] NA
This means we should no longer use proj4strings for EPSG:4326.
edzer
on 21 Sep 2020
@rsbivand, I don't believe it's related to axis order :)
@edzer, I believe I wasn't clear and the numerous breadcrumbs I left while reasoning might have been confusing, so let me summarise.
The issue is that projections that do not specify an epsg/srid input code are not supported by st_as_binary despite the fact that they might be valid and available locally.
The reason is that st_as_binary will look for the input attribute of the crs, rather than the wkt. This is what happens at L682 (the rest parses the EPSG: part if it exists.)
https://github.com/r-spatial/sf/blob/f8e6d27f3deeab543d57b1334b548dc62a6190f2/src/wkb.cpp#L680-L691
@dpprdan's original problem stems from the fact that when the file is read, the crs is set with input = "WGS 84". The wkt is fine, and it is equivalent to st_crs(4326). But because the input doesn't contain the epsg, it is lost while going through st_as_binary.
Here's a reprex with a valid existing wkt. If the crs is created with the epsg, the code goes through st_as_binary, but it is lost if it is created from wkt. I think a quick fix might be to "force" the epsg code in st_as_binary.
library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 7.0.0
x <- st_as_sfc("POINT(0 0)")
crs_epsg <- st_crs(4326)
crs_wkt <- st_crs(crs_epsg$wkt)
x %>%
st_set_crs(crs_epsg) %>%
st_as_binary(EWKB=TRUE) %>%
st_as_sfc(EWKB=TRUE) %>%
st_crs()
#> 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["unknown"],
#> AREA["World"],
#> BBOX[-90,-180,90,180]],
#> ID["EPSG",4326]]
x %>%
st_set_crs(crs_wkt) %>%
st_as_binary(EWKB=TRUE) %>%
st_as_sfc(EWKB=TRUE) %>%
st_crs()
#> Coordinate Reference System: NA
Now I understand! Thanks; this should do it.
edzer
on 21 Sep 2020
@dpprdan, thanks again for reporting this bug. I now have:
dbGetQuery(con, statement = "select st_srid(geometry) from nc_wgs84 limit 1;")
#> st_srid
#> 1 4326
with 40fb9f5.
etiennebr
on 21 Sep 2020
Thanks a lot @etiennebr and @edzer!
dpprdan
on 21 Sep 2020
Related issues
thiagoveloso
路
3Comments
Nosferican
路
3Comments
kendonB
路
4Comments
jmsigner
路
4Comments
angela-li
路
4Comments
Most helpful comment
Now I understand! Thanks; this should do it.