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

[Michael Treglia]

Wow - thanks so much for this quick fix, Edzer! I like your solution to
having syntatically different but semantically identical proj4string. Will
try this in a bit, and write back if I have any issues.
Best,
mike

On Wed, Mar 15, 2017 at 5:58 PM, Edzer Pebesma <
edzer.pebesma at uni-muenster.de> wrote:

> 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/d620f3b748d487bae910e3884c554f
> a9e3ce7090)
> 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/
>
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

	[[alternative HTML version deleted]]