forked from r-spatial/sf
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathgdal_gpkg.R
112 lines (109 loc) · 3.78 KB
/
gdal_gpkg.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
#' Set explicit timestamp to reproducibly write GeoPackage files
#'
#' Presets the timestamp for usage by GDAL by setting the environment variable
#' \code{OGR_CURRENT_DATE}.
#' After this, newly written GeoPackage files
#' created by the GDAL vector or raster driver (e.g. through
#' \code{sf::st_write()} or \code{stars::write_stars()})
#' will carry this timestamp.
#' As such the function assists in making a binary-reproducible GeoPackage file.
#'
#' The function converts the timestamp to a very specific ISO 8601 format
#' that is required by the GeoPackage standard, including conversion to UTC.
#' Cf. \href{https://www.geopackage.org/spec130/#r15}{Requirement 15} in
#' version 1.3.
#' GDAL uses the timestamp to set the \code{last_change} column of the
#' \code{gpkg_contents} table in newly written GeoPackage files.
#'
#' @param timestamp a \code{Date} or \code{POSIXct} object, used to generate
#' the timestamp.
#' For a \code{Date} object, time will be considered as \code{00:00:00} local
#' time.
#'
#' @return
#' Previous value of system variable \code{OGR_CURRENT_DATE} is returned
#' invisibly.
#'
#' @examples
#' library(openssl)
#' md5sum <- function(x) paste(md5(file(x)))
#'
#' # Using existing geopackage with vector layer:
#' filepath <- system.file("gpkg/b_pump.gpkg", package = "sf")
#' (md5_original <- md5sum(filepath))
#'
#' sf_layer <- read_sf(system.file("gpkg/b_pump.gpkg", package = "sf"))
#'
#' # A rewrite changes the checksum:
#' filepath_notimeset <- file.path(tempdir(), "b_pump_notimeset.gpkg")
#' # write 1:
#' st_write(sf_layer, dsn = filepath_notimeset)
#' (md5_notimeset1 <- md5sum(filepath_notimeset))
#' # write 2:
#' st_write(sf_layer, dsn = filepath_notimeset, delete_dsn = TRUE)
#' (md5_notimeset2 <- md5sum(filepath_notimeset))
#' # compare:
#' md5_notimeset1 == md5_notimeset2
#'
#' # Setting a fixed date
#' filepath_timeset <- file.path(tempdir(), "b_pump_timeset.gpkg")
#' (fixed_date <- as.POSIXct("2020-12-25"))
#' set_timestamp_gpkg(fixed_date)
#' # write 1 (date):
#' st_write(sf_layer, dsn = filepath_timeset)
#' md5_timeset1 <- md5sum(filepath_timeset)
#' # write 2 (date):
#' st_write(sf_layer, dsn = filepath_timeset, delete_dsn = TRUE)
#' md5_timeset2 <- md5sum(filepath_timeset)
#' # compare:
#' all.equal(md5_timeset1, md5_timeset2)
#'
#' # Setting a fixed time
#' (fixed_time <- as.POSIXct("2020-12-25 12:00:00", tz = "CET"))
#' set_timestamp_gpkg(fixed_time)
#' # write 3 (time):
#' st_write(sf_layer, dsn = filepath_timeset, delete_dsn = TRUE)
#' md5_timeset3 <- md5sum(filepath_timeset)
#' # write 4 (time):
#' st_write(sf_layer, dsn = filepath_timeset, delete_dsn = TRUE)
#' md5_timeset4 <- md5sum(filepath_timeset)
#' # compare:
#' all.equal(md5_timeset3, md5_timeset4)
#'
#' # Also works for GPKG 2D gridded coverage (with stars):
#' library(stars)
#' library(dplyr)
#'
#' filepath_stars <- file.path(tempdir(), "stars_2d.gpkg")
#' (fixed_time <- as.POSIXct("2010-02-14 12:00:00", tz = "CET"))
#' set_timestamp_gpkg(fixed_time)
#'
#' stars_2d <-
#' system.file("tif/L7_ETMs.tif", package = "stars") %>%
#' read_stars() %>%
#' slice(band, 1)
#' # write 1:
#' stars_2d %>%
#' write_stars(filepath_stars, driver = "GPKG")
#' md5_stars1 <- md5sum(filepath_stars)
#' # write 2:
#' stars_2d %>%
#' write_stars(filepath_stars, driver = "GPKG")
#' md5_stars2 <- md5sum(filepath_stars)
#' # compare:
#' all.equal(md5_stars1, md5_stars2)
#'
#' @author Floris Vanderhaeghe, \url{https://github.com/florisvdh}
#'
#' @export
set_timestamp_gpkg <- function(timestamp) {
if (!inherits(timestamp, c("Date", "POSIXct"))) {
stop("timestamp must be a Date or POSIXct object")
}
timestamp <- format(timestamp,
format = "%Y-%m-%dT%H:%M:%S.000Z",
tz = "UTC")
old <- Sys.getenv("OGR_CURRENT_DATE")
Sys.setenv(OGR_CURRENT_DATE = timestamp)
return(invisible(old))
}