JuliaClimate / ClimateTools.jl

Climate science package for Julia
https://juliaclimate.github.io/ClimateTools.jl/stable/
Other
116 stars 16 forks source link

[DOC] Improve regrid documentation #130

Closed Datseris closed 10 months ago

Datseris commented 4 years ago

Hey there,

I would like to use the regridding functionality to convert my LatLon grid data (1 degree spacing in both lon and lat) to an equal area grid. I am now writing code that produces an approximate equal area grid (which I will then immediately PR here as this seems super duper useful for me :D ).

But I find the regird function poorly documented. Firstly, what are the options for method? Secondly, may I suggest to have the regrid docstring into the "interpolation" page instead of "index"? It makes much more sense. Thirdly, this:

C = regrid(A::ClimGrid, lon::AbstractArray{N, T} where N where T, lat::AbstractArray{N, T} where N where T; dimx=[], dimy=[], method::String="linear", min=[], max=[])

is way too verbose as a docstring, and I can't imagine the information t::AbstractArray{N, T} where N where T is actually useful or relevant for the function. Fourthly, I simply don't understand the call signature regrid(A::ClimGrid, lon, lat). Which field is regrided here? The field A into the given LonLat grid?

Lastly, it should also be discussed whether grids are expected in degrees or kilometers (as some people also use kilometers).

Balinus commented 4 years ago

Hello,

thanks for the comments and suggestions. I'm gonna write more thorough information when I have some free time. Here's some quick infos:

The initial implementation was using scipy under the hood. Hence, method was refering to scipy's options (linear, bilinear and nearest neighbor I think). Now, some time ago, I quiclly re-wrote the regrid function to use GeoStats as the back-end. So, I think the method keyword is probably no-longer used. (yes, the regrid function needs love!).

The function now uses a hardcoded solution Inverse-Distance from GeoStats ecosystem: https://github.com/JuliaClimate/ClimateTools.jl/blob/b80372157c1867a6f768c064c99c28bb7c751c71/src/functions.jl#L319

My initial idea with the new system would be to provide something like

solver = Kriging(
  :precipitation => (variogram=GaussianVariogram(range=35.),)
)
R = ClimateTools.regrid(A, B, solver=solver) # regrid information from A into B coordinates

where the user provide custom solver (with a default one if not provided).

The function regrid(A::ClimGrid, lon, lat) regrid the information of ClimGrid A onto coordinates of lon and lat. I'm usually working with degrees but I'd say that if both coordinates are with the same units, GoeStats back-end should work. See this code for example (used in regrid function), where no coordinates are passed to the function. Perhaps @juliohm can confirm this assumption.

# Destination grid
SG = StructuredGrid(londest, latdest)

# Solver
n = Int(round(frac*length(lonorig[:])))
solver = InvDistWeight(target => (neighbors=n,))
SD = StructuredGridData(Dict(target => dataorig[:, :, t]), lonorig, latorig)
problem = EstimationProblem(SD, SG, target)
solution = solve(problem, solver)
OUT[:, :, t] = reshape(solution.mean[target], size(londest))
juliohm commented 4 years ago

Although the code will produce a result, (lat, lon) coordinates are special, and we definitely need more explicit support in GeoStats.jl for these coordinate systems. @Datseris is aware of this need already, and certainly this will be something to think about after we revamp the spatial data types as planned in the roadmap of the project: https://github.com/JuliaEarth/GeoStats.jl/projects

Meanwhile, we are stuck with this ugly syntax for solvers :) In the future, things should look more general, and more flexible as well. That will only happen after I finish writing 2 drafts on my TODO list for conferences this year.

consumere commented 11 months ago

Hi there, using julia 1.8 ClimateTools v0.24.1

julia> ClimateTools.regrid ERROR: UndefVarError: regrid not defined Stacktrace: [1] getproperty(x::Module, f::Symbol) @ Base ./Base.jl:31 [2] top-level scope @ REPL[184]:1

did i miss something?

Balinus commented 10 months ago

Thanks for the report, I'll take a look!

Balinus commented 10 months ago

ok, found it. I can't remember when, but regrid was renamed griddata to reflect the use of scipy.griddata under the hood. The whole approach using GeoStats was difficult the implement with enough generality on my side and there was too many spacial "use-case" from the users. So our thought was to simply remove it and let the users implement their regridding approach using GeoStats if they wanted it or to use ClimateTools.griddata for a more simpler interpolation.

I will try to update the documentation to show how GeoStats could be used.

juliohm commented 10 months ago

Thank you @Balinus for updating the docs. I would also recommend adding a link to the book in the docs for those unfamiliar with the stack we are building: https://juliaearth.github.io/geospatial-data-science-with-julia