Cropping with filter_spatial

datascience · 11 January 2022 08:36

Hi,

we have come across some questions while experimenting with the function filter_spatial. This is a reprex of the code we used:

library(openeo)
library(stars)

collections = "SENTINEL2_L2A_SENTINELHUB"
bands = c("B02", "B03", "B04", "B08", "SCL")
ext = list(11.14079 , 49.75891, 11.1432, 49.76039)
period = c("2021-04-01", "2021-06-30")
names(ext) = c("west", "south", "east", "north")
dir = tempdir()

# change to login with your credentials
# con = openeo::connect("https://openeo.cloud")
# login(
#     login_type = "oidc"
#     , provider = "egi"
#     , config = list(
#         client_id = "<client_id>"
#         , secret = "<secret>"
#     )
# )

procs = processes()

cube = procs$load_collection(
    id = collections
    , spatial_extent = ext
    , temporal_extent = period
    , bands = bands
)

# geoJSON in right format
gmj = '{"type":"Polygon","coordinates":[[[11.142508,49.759305999999998],[11.143205,49.759535],[11.142304,49.760394000000008],[11.140790999999999,49.759847],[11.141800000000002,49.758911999999998],[11.142508,49.759305999999998]]]} '
class(gmj) = "geojson"
gms = jsonlite::fromJSON(gmj)

cube_proj = procs$resample_spatial(
    data = cube
    , projection = 4326
    , method = "near"
)

cube_crop = procs$filter_spatial(
    cube_proj
    , geometries = gms
)

## create and start job
job = openeo::create_job(cube_crop)
openeo::start_job(job)

id = as.character(job$id)
jobs = openeo::list_jobs()
jobs[[id]]


## download files
openeo::download_results(
    job = jobs[[id]]$id
    , folder = dir
)

## create stars object
fls = list.files(dir, pattern = "^openEO_", full.names = TRUE)
stars = stars::read_stars(fls, along = "time")
stars = stats::setNames(stars, nm = "value")

# map
mapview::mapviewOptions(fgb = FALSE)

sf = geojsonsf::geojson_sf(gmj)
m1 = mapview::mapview(sf)
m2 = stars[,,,,6] |>  abind::adrop() |> split() |> mapview::mapview()

m1 + m2

Some things we noticed:

From the function description it seemed that the function is supposed to crop to the Polygon extent, however it looks more like it crops to the bounding box
In order for the function to work the GeoJSON string has to be manipulated to exclude the header if the polygon was part of a feature collection
What about projection? For us it worked when both the polygon and the cube were in WGS 84 (requires the cube to be converted to WGS84 first!). Is this the only option seeing as the GeoJSON does not include CRS information?

jeroen.dries · 11 January 2022 12:38

Hi,

indeed, the masking is not yet applied (using mask_polygon would work), we’ll need to add that.
This needs to be improved, we’re working on improved ‘vector-cube’ handling in any case.
Yes, unfortunately, geojson only allows for WGS84 according to the spec. So this is the safest option.

Side note: we implemented filter_spatial mainly for this use case:
https://open-eo.github.io/openeo-python-client/cookbook/sampling.html
In that case, the masking should work correctly, but for your current problem, I would go for mask_polygon, which will also be able to act as a spatial filter.

best regards,
Jeroen

datascience · 18 January 2022 12:39

Hi,

after experimenting again with filter_spatial and, as you suggested, mask_polygon we noticed some more things we thought might be helpful to you if we point them out. As always this is our reprex, very similar to above, but we added a big buffer around our study area.


library(openeo)
library(stars)

collections = "SENTINEL2_L2A_SENTINELHUB"
bands = c("B02", "B03", "B04", "B08", "SCL")
ext = list(11.00204 , 49.66901, 11.28197, 49.85029 )
period = c("2021-04-01", "2021-04-10")
names(ext) = c("west", "south", "east", "north")
dir = tempdir()

# change to login with your credentials
# con = openeo::connect("https://openeo.cloud")
# login(
#     login_type = "oidc"
#     , provider = "egi"
#     , config = list(
#         client_id = "<client_id>"
#         , secret = "<secret>"
#     )
# )

procs = processes()

cube = procs$load_collection(
    id = collections
    , spatial_extent = ext
    , temporal_extent = period
    , bands = bands
)

# geoJSON in right format
gmj = '{"type":"Polygon","coordinates":[[[11.142508,49.759305999999998],[11.143205,49.759535],[11.142304,49.760394000000008],[11.140790999999999,49.759847],[11.141800000000002,49.758911999999998],[11.142508,49.759305999999998]]]} '
class(gmj) = "geojson"
gms = jsonlite::fromJSON(gmj)

cube_proj = procs$resample_spatial(
    data = cube
    , projection = 4326
    , method = "near"
)

cube_filter = procs$filter_spatial(
    cube_proj
    , geometries = gms
)

cube_spmask = procs$mask_polygon(
    data = cube_filter
    , mask = gms
    , replacement = NA
)

# back to AOI projection

cube_reproj = procs$resample_spatial(
    cube_spmask
    , projection = 3035
    , method = "near"
)

## create and start job
job = openeo::create_job(cube_reproj)
openeo::start_job(job)

id = as.character(job$id)
jobs = openeo::list_jobs()
jobs[[id]]


## download files
openeo::download_results(
    job = jobs[[id]]$id
    , folder = dir
)

## create stars object
fls = list.files(dir, pattern = "^openEO_", full.names = TRUE)
img_crop = stars::read_stars(fls, along = "time")
img_crop = stats::setNames(img_crop, nm = "value")

# map
mapview::mapviewOptions(fgb = FALSE)

sf = geojsonsf::geojson_sf(gmj)
m1 = mapview::mapview(sf)
cube = img_crop[,,,,6] |>  abind::adrop() |> split() 
m2 =  mapview::mapview(cube)

m1 + m2

The masking seems to work well, only pixels within the area remain. But the cropping seems to only have happened on two sides of the Polygon. The rest is filled with the NAs introduced by masking. We checked that for the cube without cropping the area is located in the middle of the cube (the white dot in the middle), so this has to have happened during the cropping process.

Best!

jeroen.dries · 21 January 2022 10:32

Thanks, was this with the geotiff format?
Cropping can sometimes indeed be confusing, and for some use cases it doesn’t really matter if there’s some additional nodata, while for others this is a problem. I can certainly investigate and improve.

(Sorry for not replying earlier, I seem to have an issue with my notification mails.)

datascience · 24 January 2022 14:37

Yes exactly, we used the GeoTiff format.