Shifting coordinates with resample

[alex-s-gardner]

I have a 0.5 degree global dataset ras:

ras
360×720×192 Raster{Float64,3} with dimensions: 
  Y Projected{Float64} -89.75:0.5:89.75 ForwardOrdered Regular Points crs: EPSG,
  X Projected{Float64} -179.75:0.5:179.75 ForwardOrdered Regular Points crs: EPSG,
  Ti Sampled{DateTime} DateTime[2002-04-17T12:00:00, …, 2020-12-16T12:00:00] ForwardOrdered Irregular Points
extent: Extent(Y = (-89.75, 89.75), X = (-179.75, 179.75), Ti = (DateTime("2002-04-17T12:00:00"), DateTime("2020-12-16T12:00:00")))
missingval: NaN
crs: EPSG:4326

when I resample to 1/8 degree the longitude dimensions stay centered but the latitude dimensions sift southward

ras2 = resample(ras; res = 1/8)
1440×2880×192 Raster{Float64,3} with dimensions: 
  Y Projected{Float64} LinRange{Float64}(89.625, -90.25, 1440) ReverseOrdered Regular Intervals crs: EPSG,
  X Projected{Float64} LinRange{Float64}(-180.0, 179.875, 2880) ForwardOrdered Regular Intervals crs: EPSG,
  Ti Sampled{DateTime} DateTime[2002-04-17T12:00:00, …, 2020-12-16T12:00:00] ForwardOrdered Irregular Points
extent: Extent(Y = (-90.25, 89.75), X = (-180.0, 180.0), Ti = (DateTime("2002-04-17T12:00:00"), DateTime("2020-12-16T12:00:00")))
missingval: NaN
crs: EPSG:4326

[alex-s-gardner]

Does Rasters use pixel-is-point or pixel-is-area convention?

[rafaqz]

It depends on what the data-source specifies.

You appear to start with Points but end up with Intervals. That may not be what you want.

It also does seem like the new intervals are wrong... will need to look into it.

(If you can provide or generate a file for a full MWE that would help)

[alex-s-gardner]

So I want to initialize a 0.5 degree global grid with a time dimension:

julia> res = 1/2;

julia> lat = Y((-90+res/2):res:(90-res/2));

julia> lon = X((-180+res/2):res:(180-res/2));

julia> time = Ti(DateTime(2000,1,1):Year(1):DateTime(2020,1,1));

julia> ras = Raster(zeros(lat, lon, time), crs=EPSG(4326), missingval = NaN)
360×720×21 Raster{Float64,3} with dimensions: 
  Y Projected{Float64} -89.75:0.5:89.75 ForwardOrdered Regular Points crs: EPSG,
  X Projected{Float64} -179.75:0.5:179.75 ForwardOrdered Regular Points crs: EPSG,
  Ti Sampled{DateTime} DateTime("2000-01-01T00:00:00"):Year(1):DateTime("2020-01-01T00:00:00") ForwardOrdered Regular Points
extent: Extent(Y = (-89.75, 89.75), X = (-179.75, 179.75), Ti = (DateTime("2000-01-01T00:00:00"), DateTime("2020-01-01T00:00:00")))
missingval: NaN
crs: EPSG:4326

OK, this looks about right. I see data is being treated as points rather than cells (can't figure out how to change that) but I can work with this for now.

Now I need to increase the resolution to 1/8 degree

julia> ras2 = resample(ras; res = 1/8)
1440×2880×21 Raster{Float64,3} with dimensions: 
  Y Projected{Float64} LinRange{Float64}(89.625, -90.25, 1440) ReverseOrdered Regular Intervals crs: EPSG,
  X Projected{Float64} LinRange{Float64}(-180.0, 179.875, 2880) ForwardOrdered Regular Intervals crs: EPSG,
  Ti Sampled{DateTime} DateTime("2000-01-01T00:00:00"):Year(1):DateTime("2020-01-01T00:00:00") ForwardOrdered Regular Points
extent: Extent(Y = (-90.25, 89.75), X = (-180.0, 180.0), Ti = (DateTime("2000-01-01T00:00:00"), DateTime("2020-01-01T00:00:00")))
missingval: NaN
crs: EPSG:4326

For some reason the extents have been shifted which seems strange

[rafaqz]

Yeah there is surely a bug.

If you want intervals, pass the keywords ; sampling=Intervals(Center()) to the X and Y constructors after the ranges.

[alex-s-gardner]

I'm assuming that Intervals and Center are functions exported by DimensionalData?

[rafaqz]

No, bht they are exppted from DimensionalData.LookupArrays

[alex-s-gardner]

Looks like the issue has to be with how AG.gdalwarp is reading coordinates from the Raster as I traced the fags being passed to AG.gdalwarp and they look fine.

[rafaqz]

Ok we may need to handle it to fix anything it breaks.

[alex-s-gardner]

Another data point, if I export the original raster as a geotiff, then use GDAL at the command line to do the warping, then read back it... the extents are correct:

julia> using Rasters
julia> using Dates
julia> using ArchGDAL
julia> using DimensionalData.LookupArrays
julia> res = 1/2;
julia> lat = Y((-90+res/2):res:(90-res/2); sampling=Intervals(Center()));
julia> lon = X((-180+res/2):res:(180-res/2); sampling=Intervals(Center()));
julia> ras = Raster(zeros(lat, lon), crs=EPSG(4326), missingval = NaN)
360×720 Raster{Float64,2} with dimensions: 
  Y Projected{Float64} -89.75:0.5:89.75 ForwardOrdered Regular Intervals crs: EPSG,
  X Projected{Float64} -179.75:0.5:179.75 ForwardOrdered Regular Intervals crs: EPSG
extent: Extent(Y = (-90.0, 90.0), X = (-180.0, 180.0))
missingval: NaN
crs: EPSG:4326
julia> Rasters.write("junk.tif",ras, force=true);
shell> gdalwarp -overwrite -s_srs EPSG:4326 -t_srs EPSG:4326 -tr 0.125 0.125 -r near -of GTiff /Users/gardnera/Documents/GitHub/GRACE.jl/junk.tif /Users/gardnera/Documents/GitHub/GRACE.jl/junk_resample.tif;

Creating output file that is 2880P x 1440L.
Processing /Users/gardnera/Documents/GitHub/GRACE.jl/junk.tif [1/1] : 0Using internal nodata values (e.g. nan) for image /Users/gardnera/Documents/GitHub/GRACE.jl/junk.tif.
Copying nodata values from source /Users/gardnera/Documents/GitHub/GRACE.jl/junk.tif to destination /Users/gardnera/Documents/GitHub/GRACE.jl/junk_resample.tif;.
...10...20...30...40...50...60...70...80...90...100 - done.
julia> ras2 = Rasters.Raster("junk_resample.tif")
2880×1440 Raster{Float64,2} with dimensions: 
  X Projected{Float64} LinRange{Float64}(-180.0, 179.875, 2880) ForwardOrdered Regular Intervals crs: WellKnownText,
  Y Projected{Float64} LinRange{Float64}(89.875, -90.0, 1440) ReverseOrdered Regular Intervals crs: WellKnownText
and reference dimensions: 
  Band Categorical{Int64} 1:1 ForwardOrdered
extent: Extent(X = (-180.0, 180.0), Y = (-90.0, 90.0))
missingval: NaN
crs: 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"]]

[rafaqz]

I know what the problem is ;)

The question is when I can fix it

Either Rasters or gdal is forcing Center loci Points to be interpreted as Start loci Intervals without shifting the lookup correctly before/after

resample and warp are fairly simple and clear if you want to have a look. Somewhere or other we need to do maybeshiftlocus on the X/Y lookups to make the values the Start of the interval. If GDAL is returning Intervals, we need to switch it back to Points.

[alex-s-gardner]

I think the issue lies somewhere in AG.RasterDataset but I'm quickly getting in over my head

[felixcremer]

We are converting Points to Intervals(Start)) explicitely in line 106 of resample.Jl of the ArchGDAL extension.

[rafaqz]

Ugghh yeah that's the problem, no idea why I did that. But I think it used to run maybeshiftlocus before that so it wasn't as incorrect, and now it does it afterwards. We can just delete that whole block.

Probably just deleting this is enough:

Rasters.jl/ext/RastersArchGDALExt/resample.jl

Lines 102 to 107 in 10c6c08

	if sampling(to, XDim) isa Points
	to = set(to, dims(to, XDim) => Intervals(Start()))
	end
	if sampling(to, YDim) isa Points
	to = set(to, dims(to, YDim) => Intervals(Start()))
	end

(probably this is only tested with gdal files or something so its already a Start locus and doesn't fail)

[alex-s-gardner]

@rafaqz @felixcremer note that using Intervals(Center())) for dimensions also has the same issue so it's not unique to the Points

julia> using Rasters;

julia> using Dates;

julia> using ArchGDAL;

julia> using DimensionalData.LookupArrays;

julia> res = 1/2;

julia> lat = Y((-90+res/2):res:(90-res/2); sampling=Intervals(Center()));

julia> lon = X((-180+res/2):res:(180-res/2); sampling=Intervals(Center()));

julia> ras = Raster(zeros(lat, lon), crs=EPSG(4326), missingval = NaN);

julia> ras2 = resample(ras; res = 1/8)
1440×2880 Raster{Float64,2} with dimensions: 
  Y Projected{Float64} LinRange{Float64}(89.375, -90.5, 1440) ReverseOrdered Regular Intervals crs: EPSG,
  X Projected{Float64} LinRange{Float64}(-180.0, 179.875, 2880) ForwardOrdered Regular Intervals crs: EPSG
and reference dimensions: 
  Band Categorical{Int64} 1:1 ForwardOrdered
extent: Extent(Y = (-90.5, 89.5), X = (-180.0, 180.0))
missingval: NaN
crs: EPSG:4326
parent: ...

[alex-s-gardner]

The only flags handed to warp on line 166 of resample are:

Dict{Symbol, Any} with 3 entries:
  :t_srs => "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.017453292519…
  :tr    => [0.125, 0.125]
  :r     => :near

so this specific issue is not being caused within resample

[rafaqz]

Thanks.

We just need to thoughtfully manage the transition from Center locus to Start (for GDAL compat) for both Points and Intervals, and possibly back again too afterwards.

Something is broken somewhere in that pipeline, it just requires the time to track it down.