Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

get_pixel_line (or a helper) needs to truncate #339

Closed
mdsumner opened this issue May 11, 2024 · 4 comments
Closed

get_pixel_line (or a helper) needs to truncate #339

mdsumner opened this issue May 11, 2024 · 4 comments

Comments

@mdsumner
Copy link
Contributor

mdsumner commented May 11, 2024

A 360x180 raster has the following geotransform, in bounds -180,-90,180,90 but get_pixel_line returns values < 0 and > row/col limit.

gt <- c(-180, 1, 0, 90, 0, -1)
gdalraster::get_pixel_line(cbind(c(-190, -180, 0, 180), c(0, 100, -100, -90)), gt)

#     [,1] [,2]      ## all should be NA
#[1,]  -10   90   ## (pixel < 0, x < xmin)
#[2,]    0  -10    ## (line < 0, y > ymax) 
#[3,]  180  190  ## (line > (nrow-1), y < ymin)
# [4,]  360  180 ## (line > (nrow - 1), pixel > (ncol-1), x == xmax, y  = ymin)

I wonder if we should include a special case for x == xmax, y == ymax because they are currently out of bounds and I'm not sure that's desirable.

(Also, should they both columns be set to NA when either out of bounds ?? ... maybe not).

And, I appreciate this might be a case of having a higher level wrapper function that cleans up the output for use as row/col indices, because the logic is still mathematically sound in terms of a domain that extends past an actual dataset and that is also useful (in fact that might be the way to handle global<->local index conversion mentioned in #338).

@mdsumner mdsumner changed the title get_pixel_line needs to truncate get_pixel_line (or a helper) needs to truncate May 11, 2024
@mdsumner
Copy link
Contributor Author

The more I think about it, I don't think this is a bug and should not change.

It's a property of the geotransform, not a grid, so it's correct, and useful 🙏

@ctoney
Copy link
Collaborator

ctoney commented May 12, 2024

Right, get_pixel_line() just applies the inverse geotransform and does no bounds checking. Thanks for pointing this out. I imagined it used as a lower level function with caller doing validation of the input points first, then getting pixel/line, or at least checking the output row/cols from this function. But it wasn't clearly documented that way, and looking at it now it seems better to accept a GDALRaster object, get the geotransform from it, and do bounds checking. Addressed in #340. I went ahead and merged but I am happy to modify if you think it needs more work. Once it looks satisfactory, then I'll go back and add a test.

Behavior is unchanged for gt as numeric vector. Updated documentation at:
https://usdaforestservice.github.io/gdalraster/reference/get_pixel_line.html

@mdsumner
Copy link
Contributor Author

Excellent, I agree this is a good stance, and really good addition. Thanks!

I wonder if it should be a method on the class though, ds$get_pixel_line(xy) rather than a special-case for dataset as first arg?

I haven't done a thorough check ... standalone functions tend to be independent of instances and that seems to be a style. I don't have a strong opinion, this just occurred to me, and I guess it could exist in both forms.

@ctoney
Copy link
Collaborator

ctoney commented May 12, 2024

should be a method on the class though, ds$get_pixel_line(xy)

Yes that makes sense for it to exist in both forms. Done in #341. That PR also allows xy to be a two-column data frame that will be coerced to numeric matrix (previously took matrix input only). The documentation and examples are updated to show the class method version, and the bounds-checking output.

@ctoney ctoney closed this as completed May 22, 2024
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants