File: transform.R

package info (click to toggle)
r-cran-sf 0.9-7%2Bdfsg-5
  • links: PTS, VCS
  • area: main
  • in suites: bullseye
  • size: 6,796 kB
  • sloc: cpp: 5,333; sh: 18; makefile: 2
file content (265 lines) | stat: -rw-r--r-- 10,254 bytes parent folder | download
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
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265

chk_pol = function(x, dim = class(x)[1]) {
	PolClose = function(y) {
		if (any(head(y[[1]], 1) != tail(y[[1]], 1))) # close
			y[[1]] = rbind(y[[1]], head(y[[1]], 1))
		else if (nrow(y[[1]]) == 3) # closed, but line
			return(st_polygon(dim = dim))
		y
	}
	if (length(x) > 0 && nrow(x[[1]]) > 2)
		PolClose(x)
	else
		st_polygon(dim = dim)
}

chk_mpol = function(x) {
	cln = lapply(x, function(y) unclass(chk_pol(y, class(x)[1])))
	empty = if (length(cln))
			lengths(cln) == 0
		else
			TRUE
	# print(empty)
	st_multipolygon(cln[!empty], dim = class(x)[1])
}

sanity_check = function(x) {
    d = st_dimension(x) # flags empty geoms as NA
    if (any(d == 2, na.rm = TRUE)) { # the polygon stuff
		if (inherits(x, "sfc_POLYGON"))
        	st_sfc(lapply(x, chk_pol), crs = st_crs(x))
		else if (inherits(x, "sfc_MULTIPOLYGON"))
        	st_sfc(lapply(x, chk_mpol), crs = st_crs(x))
		else
			stop(paste("no check implemented for", class(x)[1]))
    } else
        x # nocov
}

#' Transform or convert coordinates of simple feature
#'
#' Transform or convert coordinates of simple feature
#'
#' @param x object of class sf, sfc or sfg
#' @param crs coordinate reference system: integer with the EPSG code, or character with proj4string
#' @param ... ignored
#' @param aoi area of interest, in degrees: 
#' WestLongitude, SouthLatitude, EastLongitude, NorthLatitude
#' @param pipeline character; proj4 or WKT coordinate operation, to override the default operation
#' @param reverse boolean; if \code{TRUE}, the inverse operation of the pipeline is applied
#' @param partial logical; allow for partial projection, if not all points of a geometry can be projected (corresponds to setting environment variable \code{OGR_ENABLE_PARTIAL_REPROJECTION} to \code{TRUE})
#' @param check logical; perform a sanity check on resulting polygons?
#' @details Transforms coordinates of object to new projection. Features that cannot be transformed are returned as empty geometries.
#'
#' Projecting to projections not supported by GDAL may be done by \link[lwgeom]{st_transform_proj}, part of package lwgeom.
#' @examples
#' p1 = st_point(c(7,52))
#' p2 = st_point(c(-30,20))
#' sfc = st_sfc(p1, p2, crs = 4326)
#' sfc
#' st_transform(sfc, 3857)
#' @export
st_transform = function(x, crs, ...) UseMethod("st_transform")

#' @name st_transform
#' @export
#' @examples
#' st_transform(st_sf(a=2:1, geom=sfc), "+init=epsg:3857")
#' try(st_transform(sfc, 3857, aoi = c(-280,-90,180,90)))
#' if (sf_extSoftVersion()["GDAL"] >= "3.0.0") {
#'   st_transform(sfc, pipeline =
#' 	  "+proj=pipeline +step +proj=axisswap +order=2,1") # reverse axes
#'   st_transform(sfc, pipeline =
#' 	  "+proj=pipeline +step +proj=axisswap +order=2,1", reverse = TRUE) # also reverse axes
#' }
st_transform.sfc = function(x, crs = st_crs(x), ..., 
		aoi = numeric(0), pipeline = character(0), reverse = FALSE,
		partial = TRUE, check = FALSE) {

	if (length(pipeline) == 0) {
		if (is.na(st_crs(x)))
			stop("cannot transform sfc object with missing crs")
		if (missing(crs))
			stop("argument crs cannot be missing")
	}

	crs = make_crs(crs)

	# FIXME: check for wkt here too WKT TODO:
	if (grepl("+proj=geocent", crs$proj4string) && length(x) && Dimension(x[[1]]) == "XY") # add z:
		x = st_zm(x, drop = FALSE, what = "Z")

	if (partial) {
		orig = Sys.getenv("OGR_ENABLE_PARTIAL_REPROJECTION")
		if (orig != "")
			on.exit(Sys.setenv(OGR_ENABLE_PARTIAL_REPROJECTION = orig))
		Sys.setenv(OGR_ENABLE_PARTIAL_REPROJECTION = "TRUE")
	}

	if (crs != st_crs(x)) {
		ret = structure(CPL_transform(x, crs, aoi, pipeline, reverse),
			single_type = NULL, crs = crs)
		ret = st_sfc(ret)
		if (check)
			sanity_check(ret)
		else
			ret
	} else if (length(pipeline))
		st_sfc(CPL_transform(x, crs, aoi, pipeline, reverse))
	else
		x
}

#' @name st_transform
#' @export
#' @examples
#' nc = st_read(system.file("shape/nc.shp", package="sf"))
#' st_area(nc[1,]) # area from long/lat
#' st_area(st_transform(nc[1,], 32119)) # NC state plane, m
#' st_area(st_transform(nc[1,], 2264)) # NC state plane, US foot
#' library(units)
#' set_units(st_area(st_transform(nc[1,], 2264)), m^2)
st_transform.sf = function(x, crs = st_crs(x), ...) {
	x[[ attr(x, "sf_column") ]] = st_transform(st_geometry(x), crs, ...)
	x
}

#' @name st_transform
#' @export
#' @details The \code{st_transform} method for \code{sfg} objects assumes that the CRS of the object is available as an attribute of that name.
#' @examples
#' st_transform(structure(p1, proj4string = "+init=epsg:4326"), "+init=epsg:3857")
st_transform.sfg = function(x, crs = st_crs(x), ...) {
	x = st_sfc(x, crs = attr(x, "proj4string"))
	if (missing(crs))
		stop("argument crs cannot be missing")
	crs = make_crs(crs)
	structure(st_transform(x, crs, ...)[[1]], crs = crs)
}

#' @name st_transform
#' @param type character; one of \code{have_datum_files}, \code{proj}, \code{ellps}, \code{datum}, \code{units} or \code{prime_meridians}; see Details.
#' @param path character; PROJ search path to be set
#' @export
#' @details \code{sf_proj_info} lists the available projections, ellipses, datums, units, or data search path of the PROJ library when \code{type} is equal to proj, ellps, datum, units or path; when \code{type} equals \code{have_datum_files} a boolean is returned indicating whether datum files are installed and accessible (checking for \code{conus}).
#'
#' for PROJ >= 6, \code{sf_proj_info} does not provide option \code{type = "datums"}. 
#' PROJ < 6 does not provide the option \code{type = "prime_meridians"}.
#'
#' for PROJ >= 7.1.0, the "units" query of \code{sf_proj_info} returns the \code{to_meter} 
#' variable as numeric, previous versions return a character vector containing a numeric expression.
#' @examples
#' sf_proj_info("datum")
sf_proj_info = function(type = "proj", path) {

	if (type == "have_datum_files")
		return(CPL_have_datum_files(0))

	if (type == "path")
		return(CPL_get_data_dir(FALSE))
	
	if (!missing(path) && is.character(path))
		return(invisible(CPL_set_data_dir(path)))

	if (type == "network")
		return(CPL_is_network_enabled(TRUE))

    opts <- c("proj", "ellps", "datum", "units", "prime_meridians")
    if (!(type %in% opts))
		stop("unknown type") # nocov
    t <- as.integer(match(type[1], opts) - 1)
	res = CPL_proj_info(as.integer(t))
    if (type == "proj")
		res$description <- sapply(strsplit(as.character(res$description), "\n"),
			function(x) x[1])
    data.frame(res)
}

#' @name st_transform
#' @export
st_wrap_dateline = function(x, options, quiet) UseMethod("st_wrap_dateline")

#' @name st_transform
#' @param options character; should have "WRAPDATELINE=YES" to function; another parameter that is used is "DATELINEOFFSET=10" (where 10 is the default value)
#' @param quiet logical; print options after they have been parsed?
#' @export
#' @examples
#' st_wrap_dateline(st_sfc(st_linestring(rbind(c(-179,0),c(179,0))), crs = 4326))
#' library(maps)
#' wrld <- st_as_sf(maps::map("world", fill = TRUE, plot = FALSE))
#' wrld_wrap <- st_wrap_dateline(wrld, options = c("WRAPDATELINE=YES", "DATELINEOFFSET=180"),
#'    quiet = TRUE)
#' wrld_moll <- st_transform(wrld_wrap, "+proj=moll")
#' plot(st_geometry(wrld_moll), col = "transparent")
#' @details For a discussion of using \code{options}, see \url{https://github.com/r-spatial/sf/issues/280} and \url{https://github.com/r-spatial/sf/issues/541}
st_wrap_dateline.sfc = function(x, options = "WRAPDATELINE=YES", quiet = TRUE) {
	if (is.na(st_crs(x)))
		warning("crs not set: assuming geographic coordinates")
	else
		stopifnot(st_is_longlat(x))
	stopifnot(is.character(options))
	stopifnot(is.logical(quiet) && length(quiet) == 1)
	st_sfc(CPL_wrap_dateline(x, options, quiet), crs = st_crs(x))
}

#' @name st_transform
#' @export
st_wrap_dateline.sf = function(x, options = "WRAPDATELINE=YES", quiet = TRUE) {
	st_set_geometry(x, st_sfc(CPL_wrap_dateline(st_geometry(x), options, quiet), crs = st_crs(x)))
}

#' @name st_transform
#' @export
st_wrap_dateline.sfg = function(x, options = "WRAPDATELINE=YES", quiet = TRUE) {
	st_sfc(CPL_wrap_dateline(st_geometry(x), options, quiet), crs = st_crs(x))[[1]]
}

st_to_s2 = function(x) {
	# to geocentric, spherical, unit sphere:
	st_transform(x, st_crs("+proj=geocent +a=1 +b=1 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs")) # nocov
}

#' directly transform a set of coordinates
#'
#' directly transform a set of coordinates
#' @param from character description of source CRS, or object of class \code{crs} 
#' @param to character description of target CRS, or object of class \code{crs} 
#' @param pts two-column numeric matrix, or object that can be coerced into a matrix
#' @param keep logical value controlling the handling of unprojectable points. If
#' `keep` is `TRUE`, then such points will yield `Inf` or `-Inf` in the
#' return value; otherwise an error is reported and nothing is returned.
#' @param warn logical; if \code{TRUE}, warn when non-finite values are generated
#' @param authority_compliant logical; \code{TRUE} means handle axis order authority compliant (e.g. EPSG:4326 implying x=lat, y=lon), \code{FALSE} means use visualisation order (i.e. always x=lon, y=lat)
#' @return two-column numeric matrix with transformed/converted coordinates, returning invalid values as \code{Inf}
#' @export
sf_project = function(from, to, pts, keep = FALSE, warn = TRUE, 
		authority_compliant = st_axis_order()) {

	proj_from_crs = function(x) {
		if (inherits(x, "crs")) {
			x = if (sf_extSoftVersion()["proj.4"] >= "6.0.0")
				x$wkt
			else
				x$proj4string
		}
		v = CPL_proj_is_valid(x)
		if (!v[[1]])
			stop(paste0(v[[2]], ": ", x))
		x[1]
	}
	if (!is.logical(keep) || length(keep) != 1 || is.na(keep))
		stop("'keep' must be single-length non-NA logical value")

	CPL_proj_direct(c(proj_from_crs(from), proj_from_crs(to)), as.matrix(pts), 
			keep, warn, authority_compliant)
}

#' Get or set the PROJ search paths
#' @param paths the search path to be set; omit if no paths need to be set
#' @return the search path (possibly after setting to paths)
#' @export
sf_proj_search_paths = function(paths = NULL) {
	if (!is.null(paths))
		CPL_set_proj_search_paths(as.character(paths));
	CPL_get_proj_search_paths(character(0));
}