Unexpected behaviour of DataFrame(Raster)
#778
Comments
edwardlavender
changed the title
DataFrame(Raster)DataFrame(Raster)
|
For context, I also compared the behaviour of |
|
Are you sure this is not terra::extract misbehaving? My guess is that the x and y are the boundaries and sometimes it's taking the wrong cell because of floating point error |
|
You're assuming I can read R here, I haven't used it in over 5 years. But why would DataFrame(raster) give you the same answer as terra extract? I don't get it.
|
|
Thanks for your responses! I am not sure whether this is an issue of my understanding or an issue with Rasters or terra. DataFrame(Raster) gives me a dataframe of positions and values. I assumed that when I look up the raster values at those positions, I would get the same answer. But both |
|
But in GeoArrays you're using extract. If you try calling url <- "https://raw.githubusercontent.com/edwardlavender/patter/main/inst/extdata/dat_gebco.tif"
geo <- terra::rast(url)
df <- data.frame(geo)
head(df, 50)just copy-pasting the last 5 rows here, which match z1 (not z2) in your original post vs |
|
Sounds bad. But you're going to need to make a Julia MWE so we can check it, if possible download your data in the script so it just runs |
|
But skipping Nas before you run extract seems risky, are you sure missing values are actually nas? |
|
My gut feeling is this may be because you are extracting the edges of pixels rather than centres, so floating point error is causing some to drop into the wrong pixel. We probably need to swap out LinRange to StepRangeLen (but that won't guarantee it) |
|
Tiem sorry I missed your comment. Yes that makes more sense. But I think extract should work too? |
|
Thanks for your responses! Okay that makes sense. So |
|
GDAL defines pixels by their corners, netcdf by their centres. You can change that manually with |
|
Thanks! Please can you provide a very quick example of how to use |
using Rasters
using Rasters.Lookups
# ras is some raster
# shift locus so that stored points are at the lower left of cell corners
ras2 = Rasters.shiftlocus(Start(), ras)
# shift locus so that stored points are at the center of the cell
ras3 = Rasters.shiftlocus(Center(), ras)
# shift locus so that stored points are at the top right of cell corners
ras4 = Rasters.shiftlocus(End(), ras)You can also extract "interval bounds" that are tuples of julia> Rasters.intervalbounds(cell_area, X)
1667-element Vector{Tuple{Float64, Float64}}:
(1.8479035471683413e7, 1.8479102263377886e7)
(1.8479102263377886e7, 1.8479169055072367e7)
(1.8479169055072367e7, 1.847923584676684e7)
(1.847923584676684e7, 1.8479302638461318e7)
(1.8479302638461318e7, 1.847936943015579e7)
(1.847936943015579e7, 1.847943622185027e7)
(1.847943622185027e7, 1.8479503013544742e7)
(1.8479503013544742e7, 1.847956980523922e7)
(1.847956980523922e7, 1.8479636596933696e7)
(1.8479636596933696e7, 1.8479703388628174e7)
⋮
(1.8589709309430085e7, 1.858977610112456e7)
(1.858977610112456e7, 1.8589842892819036e7)
(1.8589842892819036e7, 1.8589909684513513e7)
(1.8589909684513513e7, 1.858997647620799e7)
(1.858997647620799e7, 1.8590043267902464e7)
(1.8590043267902464e7, 1.8590110059596945e7)
(1.8590110059596945e7, 1.8590176851291418e7)
(1.8590176851291418e7, 1.8590243642985895e7)
(1.8590243642985895e7, 1.859031043468037e7)
(1.859031043468037e7, 1.8590377226374846e7)(and same for |
|
@lazarusA it would be good to add this to the docs in your PR, to show what shiftlocus means. Maybe also plot the DimPoints to really drive it home? |
|
Thanks all for the responses & example code! |
|
We should probably change (currently you can do that with
DD issue here now: |
I am confused by the behaviour of
DataFrame(Raster). For the same file, I compared the values reported byDataFrame(Raster)andR'sterra::extract(). On the same file, and for the same coordinates, there are notable discrepancies.The reprex below is written in
R, for convenient comparison ofRasterswithterra. I read the same file intoRandJulia. InJulia, I turn theRasterinto aDataFramethat I send back to R. I then compare the difference between the values on the raster reported byRastersandterra. These are up to 27 m for this example file.The text was updated successfully, but these errors were encountered: