Skip to content

Commit

Permalink
st_warp (stars native) flips longitudes a full cycle
Browse files Browse the repository at this point in the history 
when longitudes are outside the target longitude bounds, they will
be flipped a full cycle (360 degrees) in the direction of the closest bound.
st_warp(..., use_gdal=TRUE) already did this.
relevant for #256 #264 #269
  • Loading branch information
@edzer
edzer committed Mar 31, 2020
1 parent fcfa042 commit d80d98b
Show file tree
Hide file tree
Showing 2 changed files with 9 additions and 6 deletions.
5 changes: 5 additions & 0 deletions R/stars.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -198,6 +198,11 @@ colrow_from_xy = function(x, obj, NA_outside = FALSE) {
if (inherits(obj, "dimensions"))
gt = get_geotransform(obj)

if (isTRUE(st_is_longlat(st_crs(obj)))) {
bb = st_bbox(obj)
sign = ifelse(x[,1] < bb["xmin"], 1., ifelse(x[,1] > bb["xmax"], -1., 0.))
x[,1] = x[,1] + sign * 360.

This comment has been minimized.

Sign in to view

Copy link
@mmyrte

mmyrte Jan 26, 2022

Apparently, this logic also wraps around even if there is no need to - see attached image of the situation, I'm trying to st_crop a stars object using an sf object. (violet: base raster, greyscale: cropped raster using QGIS, red: polygon in the middle of the US, nowhere near the dateline)

I don't have time open a proper issue or even make a PR, sorry! Just jotting this down for the time being. For the record:

Browse[5]> bb
     xmin      ymin      xmax      ymax 
-84.98750  37.22917 -82.69583  39.03750 
Browse[5]> x
          [,1]     [,2]
[1,] -83.02522 37.87927
[2,] -82.46407 38.26894
Browse[5]> sign
[1]  0 -1
Browse[5]> x[,1] + sign * 360.
[1]  -83.02522 -442.46407

image

This comment has been minimized.

Sign in to view

Copy link
@mmyrte

mmyrte Jan 26, 2022

If anyone happens upon this for the same reason: I'm using sf::st_crop(crop=FALSE) for the time being since that avoids this call.

This comment has been minimized.

Sign in to view

Copy link
@mmyrte

mmyrte Mar 3, 2022 •
edited
Loading

To anyone finding this in the future: I tried producing a minimal reprex, but could not reproduce the condition.

This comment has been minimized.

Sign in to view

Copy link
@edzer

edzer Mar 24, 2022

Author Member

Thanks! I ran into this in #519 and modified it in e0bbe09 ; not sure it now does what it should do, and help welcome!
}
if (!any(is.na(gt))) { # have geotransform
inv_gt = gdal_inv_geotransform(gt)
if (any(is.na(inv_gt)))
Expand Down
10 changes: 4 additions & 6 deletions R/warp.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -52,7 +52,7 @@ default_target_grid = function(x, crs, cellsize = NA_real_, segments = NA) {
}

rename_xy_dimensions = function(x, dims) {
stopifnot(inherits(x, "stars"), inherits(dims, "dimensions"))
#stopifnot(inherits(x, "stars"), inherits(dims, "dimensions"))
dx = st_dimensions(x)
xxy = attr(dx, "raster")$dimensions
dimsxy = attr(dims, "raster")$dimensions
Expand All @@ -69,7 +69,7 @@ rename_xy_dimensions = function(x, dims) {
# transform grid x to dimensions target
# x is a stars object, target is a dimensions object
transform_grid_grid = function(x, target) {
stopifnot(inherits(target, "dimensions"))
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])
Expand All @@ -87,16 +87,14 @@ transform_grid_grid = function(x, target) {
d = st_dimensions(x)
# get col/row from x/y:
xy = colrow_from_xy(pts, x, NA_outside = TRUE)

from = x[[1]] #[,,1]
dims = dim(x)
index = matrix(1:prod(dims[dxy]), dims[ dxy[1] ], dims[ dxy[2] ])[xy]
index = matrix(seq_len(prod(dims[dxy])), dims[ dxy[1] ], dims[ dxy[2] ])[xy]
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,]
x[[i]] = x[[i]][index,] # FIXME: won't work for dims > 3?
dim(x[[i]]) = c(dim(target)[dxy], remaining_dims)
}
} else {
Expand Down

0 comments on commit d80d98b

Please sign in to comment.