zonal_statistics issues wrong error on CRS of data cube and input feature

[goergen95]

Hi @appelmar,

First of all, let me thank you for this really great package which despite its early stage already offers incredible useful functionality!!

I came across an error though which I currently do not understand when applying the zonal_statistics() function. I initiated both the date cube and my sf-object with EPSG:4326, however, when running the function the following error is issued:

Error: Data cube and input features have different SRSes Error in libgdalcubes_zonal_statistics(x, geom, agg_funcs, agg_bands, : Data cube and input features have different SRSes

I made a small reproducible example based on the package's sample data which produces the same error on my machine:

library(sf)
library(gdalcubes)
L8_files <- list.files(system.file("L8NY18", package = "gdalcubes"),
                       ".TIF", recursive = TRUE, full.names = TRUE)
aoi = st_read(system.file("nycd.gpkg", package = "gdalcubes"))
aoi = st_transform(aoi, st_crs("EPSG:4326"))
extent = st_bbox(aoi)
v = cube_view(srs="EPSG:4326", dy=300, dx=300, dt="P1M", 
              aggregation = "median", resampling = "bilinear",
              extent=list(left=extent[1], right=extent[3],
                          bottom=extent[2], top=extent[4], 
                          t0="2018-01-01", t1="2018-12-31"))

raster_cube(create_image_collection(L8_files, "L8_L1TP"), v) %>%
  select_bands(c("B04", "B05")) %>%
  apply_pixel("(B05-B04)/(B05+B04)", "NDVI") %>%
  zonal_statistics(aoi,
                   expr = "median(NDVI)", as_stars = TRUE)

I am not familiar with cpp thus I could not do more to troubleshoot the problem. Anyway, I also tested the above code on a docker container getting the same result. So I guess the issue is not with my local machine. Let me know if I can provide more information or testing. Darius.

[appelmar]

Thanks for the reproducible example! I get the same error and it seems that in the C++ implementation, OGRSpatialReference::IsSame() returns false. I'll have a closer look to find out why.

[appelmar]

https://github.com/appelmar/gdalcubes_R/commit/dafb64d5c9fa6771bb3e7d7ff4e46adb908cf745 now does on-the-fly transformation of input geometries if the coordinate reference system differs from the cube. However, in your example, dx and dy of v should also refer to lat lon, e.g., 0.005 instead of 300.

There is still one more issue when I try to plot the result with stars. I'll keep the issue open until this has been solved.

[goergen95]

Thanks for looking into it and yes, the x and y units are totally senseless in my example! I can confirm that the extraction itself now works seamlessly also in my use case.