-
Notifications
You must be signed in to change notification settings - Fork 95
Expand file tree
/
Copy pathwarp.R
More file actions
251 lines (238 loc) · 11.4 KB
/
warp.R
File metadata and controls
251 lines (238 loc) · 11.4 KB
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
## https://mastodon.social/@EvenRouault/111738407524055989
has_global_longitude = function(x) {
st_is_longlat(x) && isTRUE(all.equal(as.numeric(st_bbox(x)$xlim), c(-180.,180.)))
}
cut_latitude_to_3857 = function(bb) { # Pseudomercator does not have global coverage:
bb["ymin"] = max(bb["ymin"], -85.06)
bb["ymax"] = min(bb["ymax"], 85.06)
bb
}
# create a target grid from x and crs, determining cellsize if unknown
# return a dimensions object
default_target_grid = function(x, crs, cellsize = NA_real_, segments = NA) {
bb_x = st_bbox(x)
old_crs = st_crs(x)
if (st_is_longlat(x) && crs == st_crs(3857)) # common case:
bb_x = cut_latitude_to_3857(bb_x)
envelope = st_as_sfc(st_set_crs(bb_x, NA_crs_))
# global adjustment: needed to have st_segmentize span global extent
# https://github.com/r-spatial/mapview/issues/256
# https://github.com/r-tmap/tmap/issues/526 :
envelope = st_set_crs(st_segmentize(envelope, st_length(st_cast(envelope, "LINESTRING"))/segments), old_crs)
envelope_new = st_transform(envelope, crs)
bb = st_bbox(envelope_new) # in new crs
if (any(is.na(cellsize))) {
area = if (st_is_longlat(crs)) # we need a cell size in degree lon lat
diff(bb[c("xmin", "xmax")]) * diff(bb[c("ymin", "ymax")])
else
st_area(envelope_new)
ratio = if (has_rotate_or_shear(x)) {
d = st_dimensions(x)
xy = xy_from_colrow(rbind(c(0,0), c(1,1)), st_geotransform(d))
cellarea = sum((xy[1,] - xy[2,])^2) / 2 # lenght is cell diagonal
xy_dims = attr(d, "raster")$dimensions
unclass(st_area(envelope)) / (prod(dim(x)[xy_dims]) * cellarea) # > 1
} else
1.0
dxy = attr(st_dimensions(x), "raster")$dimensions
cellsize = sqrt(unclass(area)/prod(dim(x)[dxy])/ratio)
# TODO: divide by st_area(evelope_new)/st_area(envelope) ?
}
cellsize = rep(abs(cellsize), length.out = 2)
nx = ceiling(diff(bb[c("xmin", "xmax")])/cellsize[1])
ny = ceiling(diff(bb[c("ymin", "ymax")])/cellsize[2])
if (has_global_longitude(x)) { # if global coverage, don't cross the boundaries:
cellsize[1] = diff(bb[c("xmin", "xmax")])/nx
cellsize[2] = diff(bb[c("ymin", "ymax")])/ny
}
x = create_dimension(from = 1, to = nx, offset = bb["xmin"], delta = cellsize[1], refsys = crs)
y = create_dimension(from = 1, to = ny, offset = bb["ymax"], delta = -cellsize[2], refsys = crs)
create_dimensions(list(x = x, y = y), get_raster())
}
rename_xy_dimensions = function(x, dims) {
#stopifnot(inherits(x, "stars"), inherits(dims, "dimensions"))
dx = st_dimensions(x)
xxy = attr(dx, "raster")$dimensions
dimsxy = attr(dims, "raster")$dimensions
if (all(xxy == dimsxy))
x
else {
na = names(dx)
names(dx)[match(xxy, na)] = dimsxy
attr(dx, "raster")$dimensions = dimsxy
structure(x, dimensions = dx)
}
}
# transform grid x to dimensions target
# x is a stars object, target is a dimensions object
transform_grid_grid = function(x, target, threshold) {
stopifnot(inherits(x, "stars"), inherits(target, "dimensions"))
x = rename_xy_dimensions(x, target) # so we can match by name
xy_names = attr(target, "raster")$dimensions
new_pts = st_coordinates(target[xy_names])
dxy = attr(target, "raster")$dimensions
from = st_crs(target)
pts = sf::sf_project(from = from, to = st_crs(x), pts = new_pts)
# at xy (target) locations, get values from x, or put NA
# to array:
d = st_dimensions(x)
# get col/row from x/y:
xy = if (is_curvilinear(x)) {
if (!requireNamespace("FNN", quietly = TRUE))
stop("package FNN required, please install it first") #nocov
if (st_is_longlat(x))
warning("using Euclidean distance measures on geodetic coordinates")
fnn = FNN::get.knnx(cc_x <- st_coordinates(d[xy_names]), pts, 1)
if (is.na(threshold)) {
p12 = st_as_sf(as.data.frame(cc_x[1:2,]), coords = 1:2)
threshold = signif(st_distance(p12)[1,2])
message(paste("threshold set to", format(threshold),
": set a larger value if you see missing values where there shouldn't be"))
}
i = fnn$nn.index - 1
i[fnn$nn.dist > threshold] = NA
ny = dim(x)[1]
cbind(i %% ny, i %/% ny) + 1
} else
colrow_from_xy(pts, x, NA_outside = TRUE, flip = TRUE)
dims = dim(x)
index = matrix(seq_len(prod(dims[dxy])), dims[ dxy[1] ], dims[ dxy[2] ])[xy]
x = unclass(x) # avoid using [[<-.stars:
if (length(dims) > 2) {
remaining_dims = dims[setdiff(names(dims), dxy)]
newdim = c(prod(dims[dxy]), prod(remaining_dims))
for (i in seq_along(x)) {
dim(x[[i]]) = newdim
x[[i]] = x[[i]][index,]
dim(x[[i]]) = c(dim(target)[dxy], remaining_dims)
}
} else {
for (i in seq_along(x)) {
x[[i]] = x[[i]][index]
dim(x[[i]]) = dim(target)
}
}
d[dxy] = target[1:2]
st_stars(x, dimensions = create_dimensions(d, attr(target, "raster")))
}
#' Warp (resample) grids in stars objects to a new grid, possibly in an new coordinate reference system
#'
#' @param src object of class \code{stars} with source raster
#' @param dest object of class \code{stars} with target raster geometry
#' @param crs coordinate reference system for destination grid, only used when \code{dest} is missing
#' @param cellsize length 1 or 2 numeric; cellsize in target coordinate reference system units
#' @param segments (total) number of segments for segmentizing the bounding box before transforming to the new crs
#' @param use_gdal logical; if \code{TRUE}, use gdal's warp or warper, through \link[sf]{gdal_utils}
#' @param options character vector with options, passed on to gdalwarp
#' @param no_data_value value used by gdalwarp for no_data (NA) when writing to temporary file;
#' not setting this when \code{use_gdal} is \code{TRUE} leads to a warning
#' @param debug logical; if \code{TRUE}, do not remove the temporary gdalwarp destination file, and print its name
#' @param method character; see details for options; methods other than \code{near} only work when \code{use_gdal=TRUE}
#' @param threshold numeric; distance threshold for warping curvilinear grids: new cells at distances larger than threshold are assigned NA values.
#' @param ... ignored
#' @details \code{method} should be one of \code{near}, \code{bilinear}, \code{cubic}, \code{cubicspline}, \code{lanczos}, \code{average}, \code{mode}, \code{max}, \code{min}, \code{med}, \code{q1} or \code{q3}; see https://github.com/r-spatial/stars/issues/109
#' @examples
#' geomatrix = system.file("tif/geomatrix.tif", package = "stars")
#' (x = read_stars(geomatrix))
#' new_crs = st_crs('OGC:CRS84')
#' y = st_warp(x, crs = new_crs)
#' plot(st_transform(st_as_sfc(st_bbox(x)), new_crs), col = NA, border = 'red')
#' plot(st_as_sfc(y, as_points=FALSE), col = NA, border = 'green', axes = TRUE, add = TRUE)
#' image(y, add = TRUE, nbreaks = 6)
#' plot(st_as_sfc(y, as_points=TRUE), pch=3, cex=.5, col = 'blue', add = TRUE)
#' plot(st_transform(st_as_sfc(x, as_points=FALSE), new_crs), add = TRUE)
#' # warp 0-360 raster to -180-180 raster:
#' r = read_stars(system.file("nc/reduced.nc", package = "stars"))
#' r |> st_set_crs('OGC:CRS84') |> st_warp(st_as_stars(st_bbox(), dx = 2)) -> s
#' plot(r, axes = TRUE) # no CRS set, so no degree symbols in labels
#' plot(s, axes = TRUE)
#' # downsample raster (90 to 270 m)
#' r = read_stars(system.file("tif/olinda_dem_utm25s.tif", package = "stars"))
#' r270 = st_as_stars(st_bbox(r), dx = 270)
#' r270 = st_warp(r, r270)
#' @details For gridded spatial data (dimensions \code{x} and \code{y}), see figure; the existing grid is transformed into a regular grid defined by \code{dest}, possibly in a new coordinate reference system. If \code{dest} is not specified, but \code{crs} is, the procedure used to choose a target grid is similar to that of \link[raster]{projectRaster}. This entails: (i) the envelope (bounding box polygon) is transformed into the new crs, possibly after segmentation (red box); (ii) a grid is formed in this new crs, touching the transformed envelope on its East and North side, with (if cellsize is not given) a cellsize similar to the cell size of \code{src}, with an extent that at least covers \code{x}; (iii) for each cell center of this new grid, the matching grid cell of \code{x} is used; if there is no match, an \code{NA} value is used.
#' @export
st_warp = function(src, dest, ..., crs = NA_crs_, cellsize = NA_real_, segments = 100,
use_gdal = FALSE, options = character(0), no_data_value = NA_real_, debug = FALSE,
method = "near", threshold = NA_real_) {
if (!inherits(src, "stars_proxy"))
src = st_normalize(src)
if (!is.na(crs))
crs = st_crs(crs)
if (!missing(dest)) {
if (inherits(dest, "crs"))
stop("target crs should be specified with crs = ..., not as argument dest")
} else if (is.na(crs) && !use_gdal)
stop("when use_gdal is FALSE, either dest or crs need to be specified")
ret = if (use_gdal) {
if (is_curvilinear(src))
stop("for warping curvilinear grids using GDAL, use gdal_utils(\"warp\",...) directly on the source dataset of this object")
if (!is.na(no_data_value))
options = c(options, "-dstnodata", no_data_value)
else
warning("no_data_value not set: missing values will appear as zero values")
options = c(options, "-r", method)
if (all(!is.na(cellsize))) {
cellsize = rep(abs(cellsize), length.out = 2)
options = c(options, "-tr", cellsize[1], cellsize[2])
}
if (! inherits(src, "stars_proxy")) {
src = st_as_stars_proxy(src, NA_value = no_data_value)
if (debug)
cat("Writing input to: ", src[[1]], "\n")
else
on.exit(unlink(src[[1]])) # temp file
}
# collapse a multi-file proxy objects into single-file/multi-band:
if (length(src[[1]]) > 1 || length(src) > 1) {
tmp = tempfile(fileext = ".vrt") # multi-band destination
sf::gdal_utils("buildvrt", unlist(src), tmp, options = "-separate")
src[[1]] = tmp
}
if (missing(dest) && is.na(crs)) {
delete = TRUE
dest = tempfile(fileext = ".tif")
sf::gdal_utils("warp", src[[1]], dest, options = options)
} else { # dest exists, and should be used: should use warper rather than warp
# https://github.com/r-spatial/stars/issues/407
if (missing(dest)) {
dest = default_target_grid(src, crs = crs, cellsize = cellsize, segments = segments) # dimensions
dest = st_stars(list(values = array(NA_real_, dim(dest))), dest)
}
if (length(dim(src)) == 3 && length(dim(dest)) == 2)
dest = merge(do.call(c, lapply(seq_len(dim(src)[3]), function(x) dest)))
dest = if (! inherits(dest, "stars_proxy")) {
dest[[1]] = NA_real_ * dest[[1]] # blank out values
delete = !debug
st_as_stars_proxy(dest, NA_value = no_data_value)[[1]]
} else {
delete = FALSE
dest[[1]] # the file name of the stars_proxy, not to be deleted
}
sf::gdal_utils("warper", src[[1]], dest, options = method) # not "warp"!
}
if (debug)
cat("Writing result to: ", dest, "\n")
else if (delete && !inherits(src, "stars_proxy"))
on.exit(unlink(dest)) # a temp file
read_stars(dest)
} else {
if (method != "near")
stop("methods other than \"near\" are only supported if use_gdal=TRUE")
if (missing(dest))
dest = default_target_grid(src, crs = crs, cellsize = cellsize, segments = segments)
if (!inherits(dest, "stars") && !inherits(dest, "dimensions"))
stop("dest should be a stars object, or a dimensions object")
transform_grid_grid(st_as_stars(src), st_dimensions(dest), threshold)
}
# restore attributes?
if (method %in% c("near", "mode")) {
a = lapply(src, attributes)
for (i in seq_along(ret))
ret[[i]] = structure(ret[[i]],
levels = attr(ret[[i]], "levels") %||% a[[i]]$levels,
colors = attr(ret[[i]], "colors") %||% a[[i]]$colors,
class = attr(ret[[i]], "class") %||% a[[i]]$class)
}
ret
}