st_warp() with curvilinear grid part 2

[Rapsodia86]

Greetings,
I am continuing #513 issue because warping curvilinear grid gives me some problems.

So, I am exploring Sentinel 3 L2 LST (spatial res 1km) data which come in separate *.nc files.
I have followed the tutorial (https://cran.r-project.org/web/packages/stars/vignettes/stars4.html) and those steps: (#54 (comment))
to have them displayed. That is ok.
But I am getting trouble with warping to have them at the end exported to tif. See the examples later in the post.
Does anyone have any experience with those data? I am well aware of SNAP software, but I would like to avoid it as much as possible. Here is an example of processing this data using vrt in gdal (https://github.com/jgomezdans/s3_olci/blob/master/s3_olci/s3_pre_processor.py) that gives me data in tiff format (spat res: 0.009884033). I have also been exploring xarray or pyresample & kdt (based on: https://git.earthdata.nasa.gov/projects/LPDUR/repos/ecostress_swath2grid/browse/ECOSTRESS_swath2grid.py) but was never able to match SNAP output. That is why I wanted to try starts in R. Any suggestion how to properly transform to grid and export to tif using stars? Thanks!

library(stars)
> Loading required package: abind
> Loading required package: sf
> Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUE

packageVersion("stars")
>‘0.6.0’
ll_in <- ("geodetic_in.nc")
lst_in <- ("LST_in.nc")

lon <- read_stars(ll_in, sub=5) %>% adrop
lat <- read_stars(ll_in, sub=3) %>% adrop
lst <- read_stars(lst_in, sub=1) %>% adrop
ll = setNames(c(lon, lat), c("x", "y"))
temp.c = st_as_stars(lst, curvilinear = ll)
st_crs(temp.c) <- 4326

> temp.c
stars object with 2 dimensions and 1 attribute
attribute(s), summary of first 1e+05 cells:
           Min. 1st Qu.  Median     Mean 3rd Qu.   Max. NA's
LST [K] 257.794 314.068 316.592 316.0953 319.404 330.82 4101
dimension(s):
  from   to refsys                                  values x/y
x    1 1500 WGS 84 [1500x1200] 58.2526 [°],...,77.2557 [°] [x]
y    1 1200 WGS 84 [1500x1200] 28.6775 [°],...,41.9878 [°] [y]
curvilinear grid

plot(temp.c, breaks = "equal", reset = FALSE, axes = TRUE, as_points = FALSE, border = NA)

This looks ok!

#warp1 no cellsize provided
w = st_warp(temp.c, crs = 4326, no_data_value = -32768)
Warning message:
In transform_grid_grid(st_as_stars(src), st_dimensions(dest), threshold) :
  using Euclidean distance measures on geodetic coordinates
> w
stars object with 2 dimensions and 1 attribute
attribute(s), summary of first 1e+05 cells:
          Min. 1st Qu.  Median     Mean 3rd Qu.    Max.
LST [K] 247.34  304.19 315.336 313.9428 327.496 332.056
         NA's
LST [K] 11458
dimension(s):
  from   to  offset      delta refsys x/y
x    1 1580 58.2526  0.0120283 WGS 84 [x]
y    1 1140 42.3817 -0.0120283 WGS 84 [y]

Here we have an issue

#warp2 use cellsize from the previous warp1
w2 = st_warp(temp.c, crs = 4326, no_data_value = -32768, cellsize=0.012)
Warning message:
In transform_grid_grid(st_as_stars(src), st_dimensions(dest), threshold) :
  using Euclidean distance measures on geodetic coordinates
> w2
stars object with 2 dimensions and 1 attribute
attribute(s), summary of first 1e+05 cells:
           Min. 1st Qu. Median     Mean 3rd Qu.    Max.
LST [K] 294.876 317.572 325.27 322.5213 327.152 330.104
         NA's
LST [K] 97983
dimension(s):
  from   to  offset  delta refsys x/y
x    1 1584 58.2526  0.012 WGS 84 [x]
y    1 1143 42.3817 -0.012 WGS 84 [y]

That one also does not work

S3A_SL_2_LST.zip

[edzer]

Thanks, great issue! This should work better now; the key is in threshold, which now gets (better) defaults.

[Rapsodia86]

Oh great! When that commit goes thru, I will give it a try and report back:)

[edzer]

or try remotes::install_github("r-spatial/stars")

[Rapsodia86]

Ok, I just installed dev version ‘0.6.1’
I am still getting missing values:
Default approach:

w = st_warp(temp.c, crs = 4326, no_data_value = -32768)

Set cellsize based on default results (tiny differences):
w2 = st_warp(temp.c, crs = 4326, no_data_value = -32768, cellsize=0.0120283)

If using spatial res from vrt approach:
w3 = st_warp(temp.c, crs = 4326, no_data_value = -32768, cellsize=0.009884033)

And here using an image from vrt approach (spat res: 0.009884033) as a template:

#input template
temp_r
stars object with 2 dimensions and 1 attribute
attribute(s), summary of first 1e+05 cells:
             Min. 1st Qu.  Median     Mean 3rd Qu.   Max.  NA's
one5.tif  270.018 318.508 326.064 323.0789  327.66 331.11 86864
dimension(s):
  from   to  offset       delta refsys point x/y
x    1 1923 58.2529  0.00988403 WGS 84 FALSE [x]
y    1 1347 41.9878 -0.00988403 WGS 84 FALSE [y]

#output 
w4 = st_warp(temp.c,temp_r)
w4
stars object with 2 dimensions and 1 attribute
attribute(s), summary of first 1e+05 cells:
           Min. 1st Qu.  Median     Mean 3rd Qu.    Max.  NA's
LST [K] 270.018 319.618 326.276 323.5996 327.714 330.422 91951
dimension(s):
  from   to  offset       delta refsys point x/y
x    1 1923 58.2529  0.00988403 WGS 84 FALSE [x]
y    1 1347 41.9878 -0.00988403 WGS 84 FALSE [y]

[edzer]

Yes, I got it wrong, letting threshold depend on cellsize - it should depend on the spacing of temp.c. This however is irregular and expensive to compute for all neighboring pixels. Here's an attempt that uses the first two pixels only.

[Rapsodia86]

Thank you, I will check it on Monday and give a feedback!:)

[Rapsodia86]

I did an exploration on that particular file, and using the default approach (st_warp(temp.c, crs = 4326, no_data_value = -32768)) I am getting still some missing pixels with the warning (threshold set to 0.0104971091258464 : set a larger value if you see missing values where they shouldn't be). Only after I set a threshold greater than the output cellsize, do those missing values disappear. But over the upper (and lower) left corner areas, a lot of neighboring pixels have the same value. Is this expected due to the grid type?

Using another file (attached), I was setting different thresholds, but the missing values were still there.
S3B_LST.zip