library(sf)
#> Linking to GEOS 3.7.2, GDAL 2.4.2, PROJ 5.2.0
system("ogr2ogr --version", intern = TRUE)
#> [1] "GDAL 2.4.2, released 2019/06/28"
system("proj", intern = TRUE)
#> character(0)
file <- system.file("gpkg/tl.gpkg", package = "sf")
# Check proj info of original file
system(paste("ogrinfo -so -al", file), intern = TRUE)
#> [1] "INFO: Open of `/Users/ateucher/Rlibrary/sf/gpkg/tl.gpkg'"
#> [2] " using driver `GPKG' successful."
#> [3] ""
#> [4] "Layer name: tl_2016_us_state"
#> [5] "Metadata:"
#> [6] " DBF_DATE_LAST_UPDATE=2016-09-10"
#> [7] "Geometry: Polygon"
#> [8] "Feature Count: 1"
#> [9] "Extent: (-72.557100, 42.697000) - (-70.575100, 45.305800)"
#> [10] "Layer SRS WKT:"
#> [11] "GEOGCS[\"NAD83\","
#> [12] " DATUM[\"North_American_Datum_1983\","
#> [13] " SPHEROID[\"GRS 1980\",6378137,298.257222101,"
#> [14] " AUTHORITY[\"EPSG\",\"7019\"]],"
#> [15] " TOWGS84[0,0,0,0,0,0,0],"
#> [16] " AUTHORITY[\"EPSG\",\"6269\"]],"
#> [17] " PRIMEM[\"Greenwich\",0,"
#> [18] " AUTHORITY[\"EPSG\",\"8901\"]],"
#> [19] " UNIT[\"degree\",0.0174532925199433,"
#> [20] " AUTHORITY[\"EPSG\",\"9122\"]],"
#> [21] " AUTHORITY[\"EPSG\",\"4269\"]]"
#> [22] "FID Column = fid"
#> [23] "Geometry Column = geom"
#> [24] "AWATER: Integer64 (0.0)"
# convert to shp with command-line ogr2ogr
system(paste("ogr2ogr -f 'ESRI Shapefile' tl.shp", file))
# Check proj info - same as original
system("ogrinfo -so -al tl.shp", intern = TRUE)
#> [1] "INFO: Open of `tl.shp'"
#> [2] " using driver `ESRI Shapefile' successful."
#> [3] ""
#> [4] "Layer name: tl"
#> [5] "Metadata:"
#> [6] " DBF_DATE_LAST_UPDATE=2019-12-19"
#> [7] "Geometry: Polygon"
#> [8] "Feature Count: 1"
#> [9] "Extent: (-72.557124, 42.697042) - (-70.575094, 45.305778)"
#> [10] "Layer SRS WKT:"
#> [11] "GEOGCS[\"NAD83\","
#> [12] " DATUM[\"North_American_Datum_1983\","
#> [13] " SPHEROID[\"GRS 1980\",6378137,298.257222101,"
#> [14] " AUTHORITY[\"EPSG\",\"7019\"]],"
#> [15] " TOWGS84[0,0,0,0,0,0,0],"
#> [16] " AUTHORITY[\"EPSG\",\"6269\"]],"
#> [17] " PRIMEM[\"Greenwich\",0,"
#> [18] " AUTHORITY[\"EPSG\",\"8901\"]],"
#> [19] " UNIT[\"degree\",0.0174532925199433,"
#> [20] " AUTHORITY[\"EPSG\",\"9122\"]],"
#> [21] " AUTHORITY[\"EPSG\",\"4269\"]]"
#> [22] "AWATER: Integer64 (18.0)"
# read into R and check proj info - same as original
tl <- st_read(file, quiet = TRUE)
st_crs(tl)
#> Coordinate Reference System:
#> EPSG: 4269
#> proj4string: "+proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs"
# write out to new gpkg and check - it's different now
st_write(tl, "tl2.gpkg")
#> Updating layer `tl2' to data source `tl2.gpkg' using driver `GPKG'
#> Writing 1 features with 1 fields and geometry type Polygon.
system("ogrinfo -so -al tl2.gpkg", intern = TRUE)
#> [1] "INFO: Open of `tl2.gpkg'"
#> [2] " using driver `GPKG' successful."
#> [3] ""
#> [4] "Layer name: tl2"
#> [5] "Geometry: Polygon"
#> [6] "Feature Count: 1"
#> [7] "Extent: (-72.557124, 42.697042) - (-70.575094, 45.305778)"
#> [8] "Layer SRS WKT:"
#> [9] "GEOGCS[\"GRS 1980(IUGG, 1980)\","
#> [10] " DATUM[\"unknown\","
#> [11] " SPHEROID[\"GRS80\",6378137,298.257222101],"
#> [12] " TOWGS84[0,0,0,0,0,0,0]],"
#> [13] " PRIMEM[\"Greenwich\",0],"
#> [14] " UNIT[\"degree\",0.0174532925199433]]"
#> [15] "FID Column = fid"
#> [16] "Geometry Column = geom"
#> [17] "AWATER: Real (0.0)"
# Try with rgdal - sa,e thing
library(rgdal)
#> Loading required package: sp
#> rgdal: version: 1.4-8, (SVN revision 845)
#> Geospatial Data Abstraction Library extensions to R successfully loaded
#> Loaded GDAL runtime: GDAL 2.4.2, released 2019/06/28
#> Path to GDAL shared files: /Users/ateucher/Rlibrary/rgdal/gdal
#> GDAL binary built with GEOS: FALSE
#> Loaded PROJ.4 runtime: Rel. 5.2.0, September 15th, 2018, [PJ_VERSION: 520]
#> Path to PROJ.4 shared files: /Users/ateucher/Rlibrary/rgdal/proj
#> Linking to sp version: 1.3-2
tl_rgdal <- readOGR(file)
#> OGR data source with driver: GPKG
#> Source: "/Users/ateucher/Rlibrary/sf/gpkg/tl.gpkg", layer: "tl_2016_us_state"
#> with 1 features
#> It has 1 fields
#> Integer64 fields read as strings: AWATER
proj4string(tl_rgdal)
#> [1] "+proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs"
writeOGR(tl_rgdal, "tl_rgdal.gpkg", "tl", driver = "GPKG")
system("ogrinfo -so -al tl_rgdal.gpkg", intern = TRUE)
#> [1] "INFO: Open of `tl_rgdal.gpkg'"
#> [2] " using driver `GPKG' successful."
#> [3] ""
#> [4] "Layer name: tl"
#> [5] "Geometry: Polygon"
#> [6] "Feature Count: 1"
#> [7] "Extent: (-72.557124, 42.697042) - (-70.575094, 45.305778)"
#> [8] "Layer SRS WKT:"
#> [9] "GEOGCS[\"GRS 1980(IUGG, 1980)\","
#> [10] " DATUM[\"unknown\","
#> [11] " SPHEROID[\"GRS80\",6378137,298.257222101],"
#> [12] " TOWGS84[0,0,0,0,0,0,0]],"
#> [13] " PRIMEM[\"Greenwich\",0],"
#> [14] " UNIT[\"degree\",0.0174532925199433]]"
#> [15] "FID Column = fid"
#> [16] "Geometry Column = geom"
#> [17] "AWATER: String (0.0)"Created on 2019-12-19 by the reprex package (v0.3.0)