get_pixel_line (or a helper) needs to truncate

[mdsumner]

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]

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]

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]

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]

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.