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

Inconsistencies in pixel allignment after writing, aggregating and plotting heatmaps #602

Open
rafaqz opened this issue Feb 21, 2024 · 7 comments
Labels
bug Something isn't working

Comments

@rafaqz
Copy link
Owner

rafaqz commented Feb 21, 2024

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

@rafaqz rafaqz added the bug Something isn't working label May 19, 2024
@danscr
Copy link

danscr commented Sep 17, 2024

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?

@rafaqz
Copy link
Owner Author

rafaqz commented Sep 17, 2024

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
Copy link
Owner Author

rafaqz commented Sep 17, 2024

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)

@danscr
Copy link

danscr commented Sep 18, 2024

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)

@rafaqz
Copy link
Owner Author

rafaqz commented Sep 18, 2024

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
Copy link
Owner Author

rafaqz commented Sep 18, 2024

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 ┤

@danscr
Copy link

danscr commented Sep 18, 2024

Thanks for the clarification!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
bug Something isn't working
Projects
None yet
Development

No branches or pull requests

2 participants