cellsize should error with Points

[rafaqz]

Currently, Points dimensions work in cellsize and just returns very very small area values the same everywhere:

julia> using Rasters, ArchGDAL

julia> δy = 0.25
0.25

julia> lon = X(-179.875:δy:179.875);

julia> lat = Y(-89.875:δy:89.875);

julia> # Create a raster with these dimensions

julia> ras = Raster(zeros(lon, lat); crs=EPSG(4326));

julia> c_size = Rasters.cellsize(ras);

julia> extrema(c_size)
(-2.550329404864359e8, -2.550329404864359e8)

Any ideas @tiemvanderdeure ?

[tiemvanderdeure]

We should just put in a check in for this specifically. Any other issues should be caught by DimensionalData.intervalbounds

[rafaqz]

I think intervalbounds might just give the point for both sides?

Probably best to catch things here so the error is simpler for people

[tiemvanderdeure]

Yes intervalbounds just returns the point twice: https://github.com/rafaqz/DimensionalData.jl/blob/7e6d97a166ff6e5e8de2ca7e39a17f1d565dbdf3/src/Lookups/lookup_arrays.jl#L602-L603. I think that makes sense.

Apparently the floating-point math then works out so you get tiny areas everywhere.

[lazarusA]

this kinda works, but I still see negative values for the last column, so, not sure if I should trust the other outputs 😄 .

using Rasters, DimensionalData
using ArchGDAL
using IntervalSets
using DimensionalData.Lookups
δy=0.25
lat = Y(Projected(-89.875:δy:89.875;  sampling=Intervals(Start()), order = ForwardOrdered(), span = Regular(δy), crs=EPSG(4326)));
lon = X(Projected(-179.875:δy:179.875; sampling=Intervals(Start()), order = ForwardOrdered(), span = Regular(δy), crs=EPSG(4326)));

ras = Raster(zeros(lat, lon)) # lat goes first here?? what?

c_size = Rasters.cellsize(ras) # last column has negative values.

# 
c_size[X=-1..1, Y=-1..1]

[tiemvanderdeure]

You get negative values because your latitude bounds are beyond 90.

>julia intervalbounds(ras[end:end, end:end])
([(89.875, 90.125)], [(179.875, 180.125)])

You probably meant to use Intervals(Center()) here? Once you do that there are no more negative values and the total surface area matches the surface area of the earth (guess that's one very basic check).

δy=0.25
lat = Y(Projected(-89.875:δy:89.875;  sampling=Intervals(Center()), order = ForwardOrdered(), span = Regular(δy), crs=EPSG(4326)));
lon = X(Projected(-179.875:δy:179.875; sampling=Intervals(Center()), order = ForwardOrdered(), span = Regular(δy), crs=EPSG(4326)));
ras = Raster(zeros(lat, lon))
cs = cellsize(ras)
all(cs .> 0) # true
sum(cs) # very close to the earth's surface area

I don't think there are any good reasons for latitudes to be above 90 so I think this should just error. Longitudes can be beyond 180 and the sine math still works out.

lon360 = X(Projected(0.125:δy:359.875; sampling=Intervals(Center()), order = ForwardOrdered(), span = Regular(δy), crs=EPSG(4326)));
cs360 = cellsize((lat, lon360))
all(isapprox.(cs, cs360, rtol = 1e-5))

Finally, the order of the dimensions does not matter. cellsize((lat, lon)) == cellsize((lon, lat)) return true.

[lazarusA]

ohh... yes, it should be Center, I just copy/ pasted the code from the test file 😄 . This is good. Thanks a lot.