Rasters.jl icon indicating copy to clipboard operation
Rasters.jl copied to clipboard

Inconsistencies in pixel allignment after writing, aggregating and plotting heatmaps

Open rafaqz opened this issue 1 year ago • 7 comments
trafficstars

Somewhere in this pipeline there are some errors - likely with points/intervals issues and having to shift loci to the center for netcdf.

rafaqz avatar Feb 21 '24 17:02 rafaqz

I haven't gone in too much detail about the implementation of Rasters.jl, but I don't quite follow why there's an offset for aggregate. E.g., when aggregating a 100x100 raster (with X=Y=0:1:99) to scale factor 10, I'd expect a 10x10 raster with the same extent. Instead, I get a smaller extent (X' = Y' = 5:5:95). This is fine if all rasters would consider the location to be pixel center, but when writing, the GDAL geotransform considers the coordinates to be top-left corner coordinates, which results in an offset. Would this discrepancy between pixel center and corner explain the errors you mention?

danscr avatar Sep 17 '24 09:09 danscr

Aggregate can't always have the same extent, unless you have perfect integer division.

If you want to keep the size, use resample. But resample has other problems, usually you should prefer aggregate when you can.

rafaqz avatar Sep 17 '24 17:09 rafaqz

But yeah in this case seems it should work... If you have a MWE I can fix it

(Maybe you have Points? The logic may not work so well for points. An MWE will clarify this)

rafaqz avatar Sep 17 '24 17:09 rafaqz

Thanks! This should recreate the behaviour I'm seeing:

src = Raster(rand(100,100); dims=(X(0:99),Y(0:99))) #considering (0,0) as the corner
dst = aggregate(sum, src, 10) #output extent is (5, 95)

danscr avatar Sep 18 '24 07:09 danscr

Ok, what you have are Points, the default. To specify corner intervals, use this:

using Rasters.Lookups
src = Raster(rand(100,100); dims=(X(0:99; sampling=Intervals(Start()))),Y(0:99; sampling=Intervals(Start())))) #considering (0,0) as the corner
dst = aggregate(sum, src, 10)

And the bounds are ((0, 100), (0, 100))

Netcdf and GDAL have different standards here so there is no obvious default, and for now Points is it. But maybe points should become intervals after aggregation? I'm not sure.

rafaqz avatar Sep 18 '24 07:09 rafaqz

Notice how the REPL printing shows if you have points or intervals, and if intervaals have start/center/end locus.

╭───────────────────────────╮
│ 100×100 Raster{Float64,2} │
├───────────────────────────┴────────────────────────────────────────────────────────────────── dims ┐
  ↓ X Sampled{Int64} 0:99 ForwardOrdered Regular Intervals{Start},
  → Y Sampled{Int64} 0:99 ForwardOrdered Regular Intervals{Start}
├──────────────────────────────────────────────────────────────────────────────────────────── raster ┤

rafaqz avatar Sep 18 '24 07:09 rafaqz

Thanks for the clarification!

danscr avatar Sep 18 '24 08:09 danscr