resample, method=:sum errors

[lazarusA]

using DimensionalData
using DimensionalData.Lookups
using Rasters, GDAL, ArchGDAL

SINUSOIDAL_CRS = ProjString("+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs")
x_range = LinRange(-20015109.355798, 20015109.355798, 1440)
y_range = LinRange(10007554.677899, -10007554.677899, 720)

ra_data = rand(720, 1440)/1000
ra_data[rand(1:720, 350), rand(1:1440, 350)] .= NaN

raster = Raster(ra_data, (Y(y_range), X(x_range)), crs=SINUSOIDAL_CRS)
resampled = resample(raster; res=0.25, crs=EPSG(4326), method = :sum)

ERROR: GDALError (CE_Failure, code 1):
        PROJ: sinu: Invalid latitude

Stacktrace:

plus, the resampling should output an area value, isn't? due to the crs properties, which contain meter units.

[rafaqz]

For area you need:

raster = Raster(ra_data, (Y(y_range; sampling=Intervals(Start())), X(x_range; sampling=Intervals(Start())))), crs=SINUSOIDAL_CRS)

If those are in fact the start of the intervals? Probably GDAL is assuming intervals too, maybe you are using the outer bounds in the LinRange and that's why it errors?

Seems it is, this works for me:

using DimensionalData
using DimensionalData.Lookups
using Rasters, GDAL, ArchGDAL

SINUSOIDAL_CRS = ProjString("+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs")
# Here we make the ranges one step longer than we need and shorten them
# Because these are the starts of the intervals. y is 2:end because its reversed.
x_range = LinRange(-20015109.355798, 20015109.355798, 1441)[1:end-1]
y_range = LinRange(10007554.677899, -10007554.677899, 721)[2:end]

ra_data = rand(720, 1440)/1000
ra_data[rand(1:720, 350), rand(1:1440, 350)] .= NaN

raster = Raster(ra_data, (Y(y_range), X(x_range)), crs=SINUSOIDAL_CRS)
resampled = resample(raster; res=0.25, crs=EPSG(4326), method = :sum)

[rafaqz]

We should have a constructor for Raster where you pass the extent and size or resolution and it works out the intervals for you.

Like:

rast = Raster{Float64}(undef, extent; size=(X(1440), Y(721)), crs)
rast = Raster{UInt8}(undef, extent; res=(X(10.0), Y(20.0), Ti(Month(1))), crs)

[rafaqz]

Ok this is being built into create, will close this.