Rasters.jl
Rasters.jl copied to clipboard
Inconsistencies in pixel allignment after writing, aggregating and plotting heatmaps
Somewhere in this pipeline there are some errors - likely with points/intervals issues and having to shift loci to the center for netcdf.
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?
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.
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)
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)
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.
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 ┤
Thanks for the clarification!