[R-sig-Geo] Maintain SRID with st_write to postgis db

[Edzer Pebesma]

Great question! Short answer: now solved in the dev version on github;
data are written directly to postgis having epsg 2263.


Long answer:

I believe this was caused by gdal and proj.4 doing different things when
resolving epsg codes to proj4strings. sf uses proj.4 for this. When
writing a proj4string either gdal or postgis don't recognize the
proj4string that proj.4 returns on the epsg code. Neither gdal nor
postgis understand that syntactically different proj4strings may be
semantically identical.


After running your example, in postgis

# select proj4text from spatial_ref_sys where SRID = 900914;

gives:

 +proj=lcc +lat_1=41.03333333333333 +lat_2=40.66666666666666
+lat_0=40.16666666666666 +lon_0=-74 +x_0=300000 +y_0=0 +ellps=GRS80
+towgs84=0,0,0,0,0,0,0 +units=us-ft +no_defs

# select proj4text from spatial_ref_sys where SRID = 2263;

gives:

 +proj=lcc +lat_1=41.03333333333333 +lat_2=40.66666666666666
+lat_0=40.16666666666666 +lon_0=-74 +x_0=300000.0000000001 +y_0=0
+datum=NAD83 +units=us-ft +no_defs

sf causes this:

> st_crs(2263)
$epsg
[1] 2263

$proj4string
[1] "+proj=lcc +lat_1=41.03333333333333 +lat_2=40.66666666666666
+lat_0=40.16666666666666 +lon_0=-74 +x_0=300000.0000000001 +y_0=0
+ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=us-ft +no_defs"

attr(,"class")
[1] "crs"

because it uses proj.4 directly to resolve SRID strings:

/usr/share/proj/epsg has:

# NAD83 / New York Long Island (ftUS)
<2263> +proj=lcc +lat_1=41.03333333333333 +lat_2=40.66666666666666
+lat_0=40.16666666666666 +lon_0=-74 +x_0=300000.0000000001 +y_0=0
+datum=NAD83 +units=us-ft +no_defs  <>


The change I made to sf (
https://github.com/edzer/sfr/commit/d620f3b748d487bae910e3884c554fa9e3ce7090)
now first asks gdal to resolve the epsg (in case it is not NA), and
otherwise resolve the proj4string (if not NA), instead of ONLY trying to
resolve a non-NA proj4string.

On 15/03/17 20:50, Michael Treglia wrote:
> Hi All,
> 
> Been working to import and manipulate a CSV file with point data
> (lat/long), and then export to a PostGIS db.
> 
> Overall, successful, but one thing I'd like to fix - when I write out the
> layer to postgis, the SRID is not maintained. The final EPSG/SRID should be
> 2263, but when I check in PostGIS, it comes up as 900914.
> 
> Below is my code and sessionInfo, and the data are from here:
> https://data.cityofnewyork.us/Public-Safety/NYPD-Complaint-Data-Current-YTD/5uac-w243
> (downloaded as plain old CSV)
> 
> Anything I might be missing? Thanks in advance for giving a quick look!
> Mike
> 
> 
> ##Start Code
> 
> #load packages
> library(sf)
> library(RPostgreSQL)
> 
> #read data
> crime_current <- read.csv("NYPD_Complaint_Data_Current_YTD.csv",
> stringsAsFactors = FALSE)
> 
> #format time columns for easier reading in postgres (I think), as keeping
> as date format caused problems in export
> crime_current$CMPLNT_FR_TIME <-
> as.character(as.POSIXct(paste(crime_current$CMPLNT_FR_DT,
> crime_current$CMPLNT_FR_TM), format="%m/%d/%Y\ %H:%M", tz=""))
> crime_current$CMPLNT_TO_TIME <-
> as.character(as.POSIXct(paste(crime_current$CMPLNT_TO_DT,
> crime_current$CMPLNT_TO_TM), format="%m/%d/%Y\ %H:%M", tz=""))
> crime_current$RPT_DT <- as.character(as.POSIXct(crime_current$RPT_DT,
> format="%m/%d/%Y", tz=""))
> 
> #convert to sf object
> crime_current.sf <- st_as_sf(crime_current, coords = c("Longitude",
> "Latitude"), crs = 4326)
> #reproject to EPSG 2263
> crime_current.sf <- st_transform(crime_current.sf, crs=2263)
> 
> #write to postgres
> st_write(crime_current.sf, "PG:dbname=mydb user=user host=xx.xx.xx.xx",
> 'health_safety.crime_current')
> ###End Code
> 
> 
> 
>> sessionInfo()
> R version 3.3.1 (2016-06-21)
> Platform: x86_64-pc-linux-gnu (64-bit)
> Running under: Ubuntu 14.04.5 LTS
> 
> locale:
>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C
> LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8
>  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8
>  LC_PAPER=en_US.UTF-8       LC_NAME=C
>  [9] LC_ADDRESS=C               LC_TELEPHONE=C
> LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
> 
> attached base packages:
> [1] stats     graphics  grDevices utils     datasets  methods   base
> 
> other attached packages:
> [1] sp_1.2-3          RPostgreSQL_0.4-1 DBI_0.6           sf_0.3-4
> 
> loaded via a namespace (and not attached):
> [1] tools_3.3.1     units_0.4-2     Rcpp_0.12.9     udunits2_0.13
> grid_3.3.1      lattice_0.20-33
> 
> 	[[alternative HTML version deleted]]
> 
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> 

-- 
Edzer Pebesma
Institute for Geoinformatics  (ifgi),  University of Münster
Heisenbergstraße 2, 48149 Münster, Germany; +49 251 83 33081
Journal of Statistical Software:   http://www.jstatsoft.org/
Computers & Geosciences:   http://elsevier.com/locate/cageo/

-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 473 bytes
Desc: OpenPGP digital signature
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20170315/d04cdd02/attachment.bin>