rasters_jl issues

[kadyb]

Extent of the input rasters (xmax and ymax are different):

In terra:

ext(ras)
#> SpatExtent : 550785, 783915, 5611185, 5847315 (xmin, xmax, ymin, ymax)

In rasters_jl:

extent(rast)
#> Extent(X = (550785.0, 783885.0), Y = (5.611215e6, 5.847315e6), Band = (:B1, :B9))

[kadyb]

In crop.jl:

rast = Raster(RasterStack(raster_files; name=band_names, lazy=true))
x = crop(rast; to=ext)
extent(x)
#> Extent(X = (598515.0, 727485.0), Y = (5.682105e6, 5.780985e6))

xmax and ymax are different, because rasters.jl use a different option for snapping?

Extent in terra:

x = crop(ras, area, snap = "near")
ext(x)
SpatExtent : 598515, 727515, 5682105, 5781015 (xmin, xmax, ymin, ymax)

1. The cropped raster should be a deep copy in rasters.jl.

[kadyb]

In terra:

x = aggregate(ras, fact = 3, fun = "mean")
global(x[[1]], "mean", na.rm = TRUE)
#> 9767.289

x[1000,1000,1]
#> 10125.44

In rasters_jl:

x = Rasters.aggregate(mean, raster, (3))
mean(skipmissing(x[:,:,1]))
#> 2484.44657576314

x[1000,1000,1]
#> 2843.6666666666665

The number of rows and columns is also different.

[kadyb]

extract-points.jl: Everything is missing. Output should have values.

@NamedTuple{geometry::Union{Missing, Tuple{Float64, Float64}}, B1::Union{Missing, UInt16}, B10::Union{Missing, UInt16}, B11::Union{Missing, UInt16}, B2::Union{Missing, UInt16}, B3::Union{Missing, UInt16}, B4::Union{Missing, UInt16}, B5::Union{Missing, UInt16}, B6::Union{Missing, UInt16}, B7::Union{Missing, UInt16}, B9::Union{Missing, UInt16}}(((668900.0, 5.69258e6), missing, missing, missing, missing, missing, missing, missing, missing, missing, missing))

[kadyb]

In terra:

x = lapp(ras, fun = ndvi)
global(x, mean, na.rm = TRUE)
#> 0.285815

In rasters_jl:

x = get_ndvi(red, nir)
mean(skipmissing(x))
#> 0.331073887918443

[kadyb]

In write.jl NoData Value should be set to 0, not 65535.
In terra the output is saved to the temporary file.

[kadyb]

Why are the results in raster.jl different?

terra::extract(ras, buffers, fun = mean, exact = FALSE)[, -1]
#>         B1      B10      B11       B2       B3       B4       B5      B6       B7       B9
#> 1 9232.815 24792.09 23200.88 8417.715 7710.402 7387.806 11909.14 10570.8 8245.231 5051.213

terra::extract(ras, buffers, fun = mean, exact = TRUE)[, -1]
#>         B1      B10      B11       B2       B3       B4       B5       B6       B7       B9
#> 1 9234.606 24794.23 23202.35 8420.676 7714.726 7395.761 11913.44 10587.64 8257.126 5051.198

exactextractr::exact_extract(ras, buffers, fun = "mean")
#>         B1      B10      B11       B2       B3      B4       B5       B6       B7       B9
#> 1 9234.605 24794.23 23202.35 8420.676 7714.726 7395.76 11913.44 10587.64 8257.126 5051.198

zonal((Statistics.mean), rstack; of=(buffer_df.geom), progress=false)
#>  (B1 = 9242.757575757576, B10 = 24795.405594405594, B11 = 23204.291375291374, B2 = 8432.946386946387, B3 = 7730.361305361305, B4 = 7427.223776223776, B5 = 11921.039627039627, B6 = 10657.249417249417, B7 = 8306.972027972028, B9 = 5051.146853146854)

[kadyb]

Hey @asinghvi17 and @rafaqz! I ran the scripts in rasters_jl and compared the output with terra and noticed a lot of differences. Could you please check this and respond?

[asinghvi17]

Thanks for the ping. I can't run Terra on my Mac (for some reason it doesn't install) so it is great to have the results to compare to here.

For crop, what are you trying to measure? If its data loading time then let's call this something else maybe?

For the dataframe comments, julia accepts any vector of named tuples as a table, so this is completely compatible with the interface.

The differences in results are definitely something we should dig in to. I'll try to replicate those tomorrow.

[kadyb]

Thanks!

For crop, what are you trying to measure? If its data loading time then let's call this something else maybe?

I actually wrote "load", but what I mean is that after raster cropping, the output object should be a new object in memory like a deep copy (not an in-place operation). To make it comparable with other packages.

For the dataframe comments, julia accepts any vector of named tuples as a table, so this is completely compatible with the interface.

I mean that the output should be displayed as table/data frame (rows and columns).

[rafaqz]

Thanks, I will need to look through the files to see why they are different for aggregate. The extent looks like a points/intervals disagreement. nvdi may just be an incorrect implementation of the benchmark.

We can set missingval=0 in write to make that the same.

And we can add a materialisation step after crop if you want, it will at least stop that benchmark from ruining the graph as our crop is essentially free. Its designed that way to allow it's use in larger than memory work.

And yes in Julia we have the Tables.jl standard rather than a specific inbuilt dataframe object. Anything Tables.jl compatible is essentially equivalent to all other Tables jl compatible packages so it is normal to let the conversion happen at user discretion (Rasters.jl does not depend on DataFrames.jl, only Tables.jl). The Vector of NamedTuple we use is the simplest and often fastest compatible table.

In real extract workflows we don't usually use DataFrames.jl, as that's only useful if you need to use use some of its joins or filters. If you are going to write to disk it's just adding a redundant step. Running DataFrame(output) after these methods run is kind of benchmarking the DataFrames.jl package rather than Rasters.jl, but if you want to enforce DataFrame conversions that's basically all we have to add.

But it may be missing the point that the ability to use trait systems like Tables.jl are a key performance advantage of Julia over other languages. We design everything in Rasters.jl from within that context, where R packages operate in a "dataframe is core object" paradigm.

[kadyb]

Thanks! If you have any fixes, please let me know and I will run it again.

I'm unsure about the output structure for point extraction and zonal statistics tasks, because as you mentioned, in the R ecosystem it relies heavily on data frames. In the case of Python packages, I convert to pandas data frame in the scripts, so I would like the result to look the same everywhere. Can you return a table in Tables.jl that looks like data frame in R and pandas?

[rafaqz]

We can return that, but it's possible the conversion to DataFarme will double the run time!

Rasters.extract is super optimised to the point that the cost of allocating the output memory dominates the operation. Possibly DataFrames.jl can do a no-copy wrapper around our table? I will have to explore that.

[kadyb]

From what I can see, converting to data frame has no impact:

julia> @be zonal($(Statistics.mean), $rstack; of=$(buffer_df.geom), progress=false) seconds = 10
Benchmark: 76 samples with 1 evaluation
 min    122.162 ms (307961 allocs: 109.233 MiB)
 median 133.007 ms (307961 allocs: 109.233 MiB, 6.57% gc time)
 mean   132.404 ms (307961 allocs: 109.233 MiB, 5.77% gc time)
 max    136.647 ms (307961 allocs: 109.233 MiB, 8.52% gc time)

julia> @be DataFrame(zonal($(Statistics.mean), $rstack; of=$(buffer_df.geom), progress=false)) seconds = 10
Benchmark: 76 samples with 1 evaluation
 min    122.473 ms (308113 allocs: 109.386 MiB)
 median 133.389 ms (308113 allocs: 109.387 MiB, 6.77% gc time)
 mean   132.730 ms (308113 allocs: 109.387 MiB, 5.84% gc time)
 max    137.825 ms (308113 allocs: 109.387 MiB, 8.17% gc time)

But what concerns me is that the first execution takes about 5 seconds, and each subsequent one takes about 0.13 second.

[kadyb]

Setting missingval=0 in write doesn't work for me:

julia> @be write($stack_file, $ras; missingval=0, force=$true, options = $(Dict("COMPRESS" => "LZW", "INTERLEAVE" => "BAND"))) seconds=60
ERROR: GDALError (CE_Failure, code 1):
        GetNoDataValue() should be called instead

[rafaqz]

Oh good to know! I guess this is for polygons not lines or points where each row is a relatively small extract cost.

But what concerns me is that the first execution takes about 5 seconds, and each subsequent one takes about 0.13 second.

Yes, this is just the Julia language. We can get better than C performance by compiling at runtime. But at the cost that anything small that you only run once may be slow because it takes a few seconds to compile.

If your idea is to do large scale work where performance matters then compile time is irrelevant, but if you are benchmarking simple scripts you would only run once yes, the first time has some effects. So it depends on what your objective is here. We care about large scales, for example, Rasters.extract was explicitly optimised to get an hours to minutes transition for tera users, as that's the sweet spot for people who would consider using Julia for spatial work.

And with missingval your 0 is an Int64, is the file Int64? GDAL handles Int64 in a different function, and Julia treats Int as different things to floats unlike R so it's more likely to cause problems. (Either way it should be a clearer error)