Rasters.jl
Rasters.jl copied to clipboard
writing a Regular Intervals to netcdf produces Explicit Intervals
Writing a Raster to a netcdf and reading it back in again yields a Raster with Explicit Intervals, even if the raster originally had Regular Intervals.
E.g.
using Rasters, NCDatasets
x = X(-50:1:50; sampling=Rasters.Intervals());
y = Y(-50:1:50; sampling=Rasters.Intervals());
ras = Raster(rand(x, y))
write("myraster.nc", ras; force = true)
ras_nc = Raster("myraster.nc")
This isn't always trivial. For example, rasterizing a polygon to ras works, but to ras_nc does not work.
using ArchGDAL
pointvec = [
(-20.0, 30.0),
(-20.0, 10.0),
(0.0, 10.0),
(0.0, 30.0),
(-20.0, 30.0),
]
poly = ArchGDAL.createpolygon(pointvec)
rasterize(poly; to = ras, fill = true)
rasterize(poly; to = ras_nc, fill = true) # errors
Yes we should try to minimise how much we get Explicit lookups. But also we should try to convert Explicit to Regular before we rasterize as sometimes that will be possible.
Actually rasterizing into Explicit and Irregular intervals is harder and slow, because e.g. line burning is assuming a straight line accross the grid so it does the angle calculation once.
Maybe in practice most lines only cross one or a few pixels, so it wont make much difference?
I think for normal center rasterization it wont matter as much, just a level of indirection to check the current position from a vector as we step along the rows. We have very good performance so it might even still be relatively fast.
So:
- [ ] dont write
Explicitintervals if we dont have to - [x] ignore bounds matrices on read when we don't need them
- [ ] try to convert
Explicitto regular in rasterize/mask/zonal etc where this matters - [ ] write rasterization and line-burning algs that can handle irregular grids
Not entirely sure if this is totally fixed, but the MWE I gave no longer errors
I think getting Regular instead of Explicit on load will fix 99% of the problem