Comparing R Terra and Julia Rasters crop function: alignment y to the geometry of x

[Rapsodia86]

Hello!
I have my code written in R using terra package (https://github.com/rspatial/terra), where I get pixel count based on different masks from rasters & polygons.

In a simple experiment, I cropped a raster with a smaller one. (Pixels are the same resolution and they 100% overlap). And so, due to differences in crop approach, I am getting different results.

In terra's crop function, as a default, there is a parameter snap set as "near" for the alignment y to the geometry of x. I do not care about touches because it is for polygons only. (https://rdrr.io/cran/terra/man/crop.html)
Depending on the settings, I get:

(1) snap="in"
dimensions  : 901, 1760, 1  (nrow, ncol, nlyr)
extent      : -1986126, -1106126, 1312602, 1763102  (xmin, xmax, ymin, ymax)
(2) snap="out"
dimensions  : 903, 1760, 1  (nrow, ncol, nlyr)
extent      : -1986126, -1106126, 1312102, 1763602  (xmin, xmax, ymin, ymax)
(3) snap="near" #This is a default approach I use in my original code, and this I what I would like to have in Julia.
dimensions  : 902, 1760, 1  (nrow, ncol, nlyr)
extent      : -1986126, -1106126, 1312102, 1763102  (xmin, xmax, ymin, ymax)

In Julia Rasters, the control is via touches:
If touches=false:

1760×901 Raster{Float32,2} with dimensions: 
  X Projected{Float64} LinRange{Float64}(-1.98613e6, -1.10663e6, 1760) ForwardOrdered Regular Intervals crs: WellKnownText,
  Y Projected{Float64} LinRange{Float64}(1.7626e6, 1.3126e6, 901) ReverseOrdered Regular Intervals crs: WellKnownText     
and reference dimensions: 
  Band Categorical{String} String["TSen_Slope"] Unordered
extent: Extent(X = (-1.986125753830057e6, -1.1061257538300573e6), Y = (1.3126021807256853e6, 1.7631021807256853e6))       
missingval: NaN32

If touches=true:

1762×903 Raster{Float32,2} with dimensions: 
  X Projected{Float64} LinRange{Float64}(-1.98663e6, -1.10613e6, 1762) ForwardOrdered Regular Intervals crs: WellKnownText,
  Y Projected{Float64} LinRange{Float64}(1.7631e6, 1.3121e6, 903) ReverseOrdered Regular Intervals crs: WellKnownText     
and reference dimensions: 
  Band Categorical{String} String["TSen_Slope"] Unordered
extent: Extent(X = (-1.9866257538300573e6, -1.105625753830057e6), Y = (1.3121021807256853e6, 1.7636021807256856e6))       
missingval: NaN32

So the comparable result (based on pixel dimension) is when in terra I set snap="in" while in Julia Raster's touches=false. But all my original R scripts are based on default snap="near".
Also, I have seen #312 for boundary using mask.
Would you consider adding alignment control?

Thank you!

[rafaqz]

It seems terra does not move the 1760 at all?

Essentially we will never recreate exactly what terra does for all commands in Rasters.jl. stars in R doesnt have snap for crop either, and rasterio in python has all_touched just like we do.

Maybe explain a little more why using near matters to you here.

(The touches argument is to make sure boundary=touches in rasterize/mask would not have any pixels cut off. So yes, usually slighlty larger)

[Rapsodia86]

"It seems terra does not move the 1760 at all?"
Hmm, did you mean 1762 perhaps as in Rasters.jl crop with touches:true?

"Maybe explain a little more why using near matters to you here."
For a very practical reason, which is consistency. For the previous research project, I used R, and I would like to use the same approach in my new investigation. If I now change the way of extracting and then counting pixels, I will introduce some noise and bias, and that may impact my results and comparison. Besides, that is the default option which probably is not very frequently changed when someone uses two rasters to crop.

Moreover, using Julia and your package I can get results in minutes instead of hours when using R terra. That is why I wanted to do my research in Julia using Rasters. But, as it seems crop(), here has a different approach so I will have to stick to R in that part at least. And that is ok.

In terms of touches, at least in terra, it looks more like your boundary in e.g. extract() "boundary: for polygons, include pixels where the :center is inside the polygon, where the polygon :touches the pixel, or that are completely :inside the polygon. The default is :center"

"Essentially we will never recreate exactly what terra does for all commands in Rasters.jl. stars in R doesnt have snap for crop either, and rasterio in python has all_touched just like we do."

Ok, I think that explains everything and closes the topic:)
Thanks!

[rafaqz]

I think you are right, it actually makes sense to have the :center option here too, which is essentially what you want and should be roughly equal to :near. That would make this more consistent with other methods. I'm happy to make this change.

But I need to emphasize that there is no guarantee you will always get the same result as terra even after that. Points that are on the line are treated differently per axis in gdal/R but consistently here. So results are not always identical.

[Rapsodia86]

That is great! If you need help in e.g. testing I am all in!

"But I need to emphasize that there is no guarantee you will always get the same result as terra even after that. Points that are on the line are treated differently per axis in gdal/R but consistently here. So results are not always identical."

Yes, I know. And that is a general precaution, not only for terra. That is why I am testing each step and exploring!
Thank you very much!

[Rapsodia86]

Hello,
Any updates on this maybe?:)

[rafaqz]

Ahh sorry no. I suggest you make a new issue that writes out exactly what the change should be with every detail specified.

With my workload Im unfirtunately unlikely to read back over a long conversation and I just skim the issue titles as a reminder of what Im meant to do.

[rafaqz]

Closed in favor of #593