Flexible interpolation ("regridding") using YAC #546

ckhroulev · 2024-04-25T23:27:09Z

This PR contains changes needed to implement flexible interpolation. The goal is to reduce the amount of pre-processing needed to run PISM with some publicly available datasets, e.g. bed topography and RCM outputs.

Once this work is finished, users will be able to

Conservatively regrid (-regrid_file XXX.nc -regrid_vars YYY,ZZZ) from a file that uses a fine grid (e.g. a BedMachine topography dataset) to a coarser grid used by a simulation. Here conservative interpolation is needed to avoid grid aliasing issues.
Conservatively regrid climate forcings. Here conservative interpolation is needed to preserve the amount of mass added to the system due to a forcing.
Regrid from a file using a different projection.

The interpolation part of this PR is more-or-less complete. I may need to take care of some corner cases, though.

Use PISM's old interpolation code if the source and target grids match: no need to use fancy interpolation on a sphere if no interpolation is needed.
Consider supporting global longitude,latitude grids. These are different in at least two ways: a) grid cell boundaries should follow lines of constant longitude and latitude instead of great circles and b) they are cyclic, which requires special treatment.

Most of the remaining work is related to getting projection information from an input file, including the input file used as a basis for PISM's internal 2D grid.

I decided to support input files containing CF-style "grid mapping" variables and added code processing these.

Checklist

update documentation
update CHANGES.rst
add regression tests

…tput argument

It reads valid range information from a file, which is then used to check if data in this same file falls within this valid range. I don't know why I wrote this code: right now I think we should warn about valid ranges *we* know about, e.g. minimum of 0 for ice thickness -- and not others.

…orting

This will make it possible to allocate InputInterpolationXXX from within Grid without messing with shared pointers and worrying about ownership.

…ctly

…as to be a CRS See https://stackoverflow.com/questions/71823105/migrating-to-proj-api-6-without-forcing-crs-input A better option might be to try to avoid using proj_create_crs_to_crs().

This is needed so that PISM's C++ code does not have include YAXT's header files, which are C-only.

We cannot create more than one grid called "internal".

…time

…ing" variable

…f a global attribute

Also: better error messages.

…apping variable

bueler · 2024-04-26T03:04:33Z

Amazing! 🎉

ckhroulev added 30 commits April 25, 2024 13:38

Put YACInterpolation into PISM's source tree

6c4f82a

Stop if projection of the internal grid is not known

e4107af

Whitespace changes

a761437

Remove mapping(file, unit_system) from the header

21bd023

Move YACInterpolation into the "pism" namespace

c440e17

Move LonLatCalculator into YACInterpolation

c33afa0

Make sure YACInterpolation is built only when YAC and PROJ are available

5e1cc09

Clean up regridding code a bit more

0cc8446

Fix a typo in a comment

bba2f9f

The check_for_missing_values() call is redundant here

092047d

Add a documenting comment

3e4d0b1

Whitespace changes

6a4c3a0

io::read_attributes() returns VariableMetadata instead of using an ou…

8e61f12

…tput argument

Fix compiler warnings in YACInterpolation.cc

c197f06

Rename File::find_{dimension,variable} to {dimension,variable}_exists

dc9deb9

Improve VariableMetadata::report_range()

278870f

Cleaning up the interpolation code

cc378c4

Separate reading/interpolation, input validation, and input range rep…

ca530bb

…orting

Making YACInterpolation interface easier to use...

abc8b3c

Add Interpolation2D, a base class for "regridding" methods

817aa91

io::regrid_spatial_variable(): use a const reference argument

e6d47de

Add Interpolation2DRegular

f3a6a2a

Start using old interpolation code wrapped in a new API

fcf1b86

Rename pism::Interpolation -> Interpolation1D

087a314

Re-organize new input interpolation code

41bb018

Add the ability to read a specific record from an input file

06569ae

Use new interpolation code in array::Forcing

6dfbd59

Fix some comments

853ae58

Rename YACInterpolation for consistency

5dc5ad6

ckhroulev added 28 commits April 25, 2024 13:38

Don't use shared_ptr<const Grid> in InputInterpolationXXX

55b5174

This will make it possible to allocate InputInterpolationXXX from within Grid without messing with shared pointers and worrying about ownership.

Add Grid::get_interpolation()

96355ae

Use Grid::get_interpolation()

f78567d

Try to get projection info from both "proj" and "proj4" attributes

4a550f6

More descriptive messages about interpolation directions

d87b51d

Improve stdout reporting (add a space to line things up)

7393bc0

Use get_projection_info() instead of reading the "proj" argument dire…

bed54ec

…ctly

Remove "+init" from PROJ strings "+init=epsg:XXXX"

6a9257f

Get projection info from the input file *before* regridding

bba38f0

Add logging messages before potentially time consuming computations

934b216

Improve logging messages

43e9814

Store initialized 2D interpolation methods to use later

d452fd9

Slightly improve two error messages

1bf9e59

Try to work around the fact that with new(-er) PROJ the proj string h…

9fdff4e

…as to be a CRS See https://stackoverflow.com/questions/71823105/migrating-to-proj-api-6-without-forcing-crs-input A better option might be to try to avoid using proj_create_crs_to_crs().

Use strtol() since pism::parse_integer() is too restrictive

5b2ab11

Add +type=crs to strings that look like PROJ.4 strings

7856dcc

Add a YAXT "wrapper"

a812931

This is needed so that PISM's C++ code does not have include YAXT's header files, which are C-only.

Target grid names have to be unique

9744308

We cannot create more than one grid called "internal".

Initialize YAXT so we can use more than one YAC instance at the same …

35d9c58

…time

Add conversions from CF "mapping" variables to PROJ strings

87a8807

Set the "grid_mapping" attribute if an output file contains the "mapp…

0fc43a7

…ing" variable

Work on getting the grid mapping from a particular variable instead o…

d8fd47c

…f a global attribute

Get grid mapping info from the variable in the "grid_mapping" attribute

dd077f3

Move grid_name() from InputInterpolation to the pism namespace

2612a81

Don't use YAC if PISM's grid projection is not known

8ee03ae

Also: better error messages.

Report which interpolation method is used to read data from a file

d20aebc

get_projection_info() uses the name of a 2D/3D variable, not the CF m…

4fb802e

…apping variable

Improving get_projection_info()

198e0ed

ckhroulev self-assigned this Apr 25, 2024

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Flexible interpolation ("regridding") using YAC #546

Flexible interpolation ("regridding") using YAC #546

ckhroulev commented Apr 25, 2024

bueler commented Apr 26, 2024

Flexible interpolation ("regridding") using YAC #546

Are you sure you want to change the base?

Flexible interpolation ("regridding") using YAC #546

Conversation

ckhroulev commented Apr 25, 2024

bueler commented Apr 26, 2024