-
Notifications
You must be signed in to change notification settings - Fork 37
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Shifting coordinates with resample #538
Comments
Does Rasters use pixel-is-point or pixel-is-area convention? |
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) |
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 |
Yeah there is surely a bug. If you want intervals, pass the keywords |
I'm assuming that |
No, bht they are exppted from DimensionalData.LookupArrays |
Looks like the issue has to be with how |
Ok we may need to handle it to fix anything it breaks. |
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"]] |
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
|
I think the issue lies somewhere in AG.RasterDataset but I'm quickly getting in over my head |
We are converting Points to Intervals(Start)) explicitely in line 106 of resample.Jl of the ArchGDAL extension. |
Ugghh yeah that's the problem, no idea why I did that. But I think it used to run Probably just deleting this is enough: Rasters.jl/ext/RastersArchGDALExt/resample.jl Lines 102 to 107 in 10c6c08
(probably this is only tested with gdal files or something so its already a Start locus and doesn't fail) |
@rafaqz @felixcremer note that using 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: ... |
The only flags handed to warp on line 166 of 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 |
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. |
I have a 0.5 degree global dataset
ras
:when I resample to 1/8 degree the longitude dimensions stay centered but the latitude dimensions sift southward
The text was updated successfully, but these errors were encountered: