R/warp.R

In r-spatial/stars: Spatiotemporal Arrays, Raster and Vector Data Cubes

## 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)
	if (st_is_longlat(x) && crs == st_crs(3857)) # common case:
		bb_x = cut_latitude_to_3857(bb_x)
	envelope = st_as_sfc(bb_x)
	# global adjustment: needed to have st_segmentize span global extent
	# https://github.com/r-spatial/mapview/issues/256
	envelope = if (!is.na(segments) && !has_global_longitude(x)) # FIXME: should this branch be retained?
				st_segmentize(envelope, st_length(st_cast(envelope, "LINESTRING"))/segments)
			else {
				# https://github.com/r-tmap/tmap/issues/526 :
				old_crs = st_crs(envelope)
				st_crs(envelope) = NA_crs_
				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)
	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)
			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
}