Closed rafaqz closed 1 month ago
We should just put in a check in for this specifically. Any other issues should be caught by DimensionalData.intervalbounds
I think intervalbounds
might just give the point for both sides?
Probably best to catch things here so the error is simpler for people
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.
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]
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
.
ohh... yes, it should be Center
, I just copy/ pasted the code from the test file 😄 . This is good. Thanks a lot.
Currently,
Points
dimensions work incellsize
and just returns very very small area values the same everywhere:Any ideas @tiemvanderdeure ?