Subsetting causes small errors in geographic extent

[biskwikman]

I'm having an issue where when I try to subset a global dataset, the output is slightly different than expected.
Right now I'm using an ENVI file of process MODIS data with 12 bands, 1 for each month.

# My file. There is also a header file associated with it which has the geographic information.
filename = "MOD11A2.061.LST_Day.GLOBAL.30km.2021.degC.mon.bsq.flt"

A = Raster(filename)

The output of A is:

1440×720×12 Raster{Float32,3} with dimensions: 
  X Projected{Float64} LinRange{Float64}(-180.0, 179.75, 1440) ForwardOrdered Regular Intervals crs: WellKnownText,
  Y Projected{Float64} LinRange{Float64}(89.75, -90.0, 720) ReverseOrdered Regular Intervals crs: WellKnownText,
  Band Categorical{Int64} 1:12 ForwardOrdered
extent: Extent(X = (-180.0, 180.0), Y = (-90.0, 90.0), Band = (1, 12))missingval: -9999.0f0crs: GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]
parent:
[:, :, 1]
             89.75     89.5     89.25     89.0  …    -89.5       -89.75      -90.0
 -180.0   -9999.0   -9999.0  -9999.0   -9999.0     -9999.0     -9999.0     -9999.0
    ⋮                                           ⋱                          
  179.5   -9999.0   -9999.0  -9999.0   -9999.0       -31.6266    -31.1014  -9999.0
  179.75  -9999.0   -9999.0  -9999.0   -9999.0       -31.5882  -9999.0     -9999.0
[and 11 more slices...]

As you can see, it is a series of 12 1440 x 720px images.

If I try to subset it like so:

A_cropped = A[X(60 .. 180), Y(-10 .. 80)]

I believe mathmatically, as each pixel is 0.25 degrees, the result should be a series of 480 x 360. However, the resulting sizes are 479 x 360:

479×360×12 Raster{Float32,3} with dimensions: 
  X Projected{Float64} LinRange{Float64}(60.25, 179.75, 479) ForwardOrdered Regular Intervals crs: WellKnownText,
  Y Projected{Float64} LinRange{Float64}(79.75, -10.0, 360) ReverseOrdered Regular Intervals crs: WellKnownText,
  Band Categorical{Int64} 1:12 ForwardOrdered
extent: Extent(X = (60.24999999999999, 180.0), Y = (-9.999999999999993, 80.0), Band = (1, 12))missingval: -9999.0f0crs: GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]
parent:
[:, :, 1]
            79.75     79.5     79.25  …     -9.25     -9.5     -9.75    -10.0
  60.25  -9999.0   -9999.0  -9999.0      -9999.0   -9999.0  -9999.0   -9999.0
   ⋮                                  ⋱      ⋮                        
 179.5   -9999.0   -9999.0  -9999.0      -9999.0   -9999.0  -9999.0   -9999.0
 179.75  -9999.0   -9999.0  -9999.0      -9999.0   -9999.0  -9999.0   -9999.0
[and 11 more slices...]

The Extent is also listed at 60.2499999999999 by 180.0, where it should just be 60 by 180.

Here is the file and associated header (required to be in the same directory when loaded) I used.
google drive link

[rafaqz]

It cant cut pixels in half though right?

You are selecting an integer interval accross intervals that are offset by 0.25. But it can only return whole intervals (e.g. pixels)

Oh each pixel is 0.25 I get it.

Maybe your data is points or you want it to act like points?

But then you will get a different kind of error.

Ok the problem may be the LinRange floating point error with something like that .99999 I can see. Probably your axis needs to be a StepRange. I guess the LinRange is my fault from somewhere or other.

[rafaqz]

Maybe try to set it as a StepRange with e.g

A = set(A, Y => 60.0:0.25:180)

Or whatever the numbers should be.

StepRange has a TwicePrecision hack that fixes these fp errors, but also makes it kinda slow. But it should be used everywhere so this doesnt happen.

(And sorry the phone app is making a mess of my messages)

[rafaqz]

It turns out this is a hard problem... because using a StepRange occasionally causes the opposite but much more fatal error - the dimension comes out with the wrong length. More fatal because the raster wont even construct.

We could possibly catch length errors and switch to a LinRange for those cases. But unfortunately its not just a matter of swapping things out, and I'm not sure about the best design for this. Floating point error is an unfortunate fact of life.