Skip to content
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

Remap ICON to swiss #44

Merged
merged 13 commits into from
Oct 23, 2024
Merged

Remap ICON to swiss #44

merged 13 commits into from
Oct 23, 2024

Conversation

cfkanesan
Copy link
Collaborator

@cfkanesan cfkanesan commented Sep 6, 2024

Purpose

Provide a single step interpolation operator from ICON native grid to regular grids in the Swiss coordinate systems. Note that technically any regular grid in a projected, Cartesian coordinate system is admissible but warning is emitted if the coordinate system is not Swiss. The interpolation is done on a UTM that is expected to be usable over the entire ICON domain and thus any relevant output CRS can be used (rotlatlon, geolatlon, swiss, etc...)

Note that the unit is disabled due to the fieldextra implementation taking too much time to execute.

Code changes:

  • Added regrid.iconremap operator.

Requirements changes:

  • Added pyproj and scipy as dependencies.

return indices, weights


def _cropped_domain(xy, uv, buffer=1000):
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think it would make sense to provide types and a bit more doc on these functions.
Even if declare private, for maintainability and increase its readability by developers.
Arguments are passed to other internal functions and complexity is not so trivial anymore.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Also the function is not just cropping the domain, since it is also computing the weights?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, it's essentially cropping the input domain before passing it to the computation of the weights to avoid creating triangles where no interpolation is done.

return idx[indices], weights


def iconremap(
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

maybe a comment on the API:
we have

  • icon2geolatlon: implicit method (RBF), dst grid can not be passed
  • icon2rotlatlon: implicit method (RBF), dst grid can not be passed
  • iconremap: method passed as argument (byc), dst grid passed by argument

could we homogenise the api, maybe also 1 single entry point that dispatches according to the method?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't think that multiplexing the function through a single interface is a good design. Especially in this case where the behavior change between the methods is so large. On the one hand, the weights are pre computed so the destination grid is dictated by that and on the other the weights are computed on the fly so a destination grid can be defined.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If you are not considering refactoring the api, I would at least make it explicit that pre-computed weights from RBF are used, either through docstrings or by renaming those methods.

def iconremap(
field: xr.DataArray, dst: RegularGrid, method: Literal["byc"] = "byc"
) -> xr.DataArray:
"""Remap ICON native grid data to the swiss grid.
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is it necessarily the swiss grid?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

No it isn't. I'll change this

dst : RegularGrid
A regular grid in the swiss coordinate system.
method : Literal["byc"]
Method used to perform the interpolation.
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Missing description of the interpolation method.

xy = np.array(points_src).T
uv = np.array(points_dst).T

if method == "byc":
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I would simply put

if method not in ["byc"]:
    raise NotImplementedError(f"method: {method} is not implemented")

at the start of the function.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I was hoping to add more than one but that would be for another PR

Copy link
Collaborator

@frazane frazane left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

lgtm

Copy link
Contributor

@cosunae cosunae left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

looks good, thanks

@cfkanesan cfkanesan merged commit 2de56d3 into main Oct 23, 2024
2 checks passed
@cfkanesan cfkanesan deleted the regrid-linear branch October 23, 2024 13:12
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants