remapping weights script

[anton-seaice]

Contributes to #44. The script makes remapping weights which map from all land cells into the nearest ocean cell and conserves the total volume. There is no spreading of runoff, or effort to put the runoff at the mouth of a bay/river when it is further away than the closest ocean.

Note:

There is not enough memory on the gadi login node to run this, its simplest to run in a terminal through are.nci.org.au

I get this warning, but it works ok:

WARNING: There was an error initializing an OpenFabrics device.

  Local host:   gadi-cpu-clx-1209
  Local device: mlx5_0

[aekiss]

Here, does "nearest" mean "smallest difference in grid index"? This may not correspond to smallest distance on the sphere, especially at the join between the eastern and western edges, and in the tripolar region.

[anton-seaice]

centerCoords is the geographic coordinate of the center of the cell. This function returns the index of the cell with the smallest difference in the geographic coordinate (by euclidean distance, not surface of a sphere). It won't wrap from east to western edge but should be ok for everywhere else.

[aekiss]

ok, might need to take a look at whether anything fishy happens at the join - I hope there are no major rivers near 80E in the Indian or Arctic ...

[dougiesquire]

Using the haversine distance (with a small bug fix) successfully maps across the join - see #45 (review)

[aekiss]

Thanks @anton-seaice, looks good. My only significant concern is whether this is nearest-neighbour on the sphere or on the grid.

[anton-seaice]

Thanks @anton-seaice, looks good. My only significant concern is whether this is nearest-neighbour on the sphere or on the grid.

Its nearest-neighbour by euclidean distance between the center coordinates (in degrees) of the grid. This means that for the non-square cells, the results would be a bit weird (e.g. a nominal 0.25deg cell might be 6km in x and 11km in y in the Arctic, but be around 0.25 degrees in both directions).

[anton-seaice]

@aekiss

It looks like using a Ball Tree instead of a KDTree allows for 'haversine' distances

See https://scikit-learn.org/stable/modules/neighbors.html#ball-tree

and
https://scikit-learn.org/stable/modules/generated/sklearn.metrics.pairwise.haversine_distances.html#sklearn.metrics.pairwise.haversine_distances

The haversine distances are appropriate for surface of the sphere, and should wrap east to west around the globe. Its a fairly simple change by the looks of it.

[aekiss]

Cool, I don't think the difference between nearest-by-degrees and by-distance is likely to make much difference, but we may as well try use ball tree if it's easy.

[aekiss]

Oh, I just re-read - if it also maps at the east-west join, then we should definitely use it. Will it also handle mapping across the tripole seam?

[anton-seaice]

I made the update - the new scheme is slower to make weights but otherwise very similar (and still conserves runoff in the model)

[dougiesquire]

Thanks @anton-seaice!

I created a plot to convince myself that the BallTree query is working as expected and it revealed a small bug.

In the following, crosses shows every 50th cell centre. Blue crosses are land cells, yellow crosses are ocean cells. Red lines show which land cells are mapped to which ocean cells. Clearly there are land cells in the Antarctic being mapped to ocean cells in the Arctic.

As we worked out in our chat, this is because the coordinates need to be provided to the BallTree routines in the order (lat, lon), not (lon, lat). Switching this, as suggested in my review, fixes the issue.

[anton-seaice]

Thanks @dougiesquire - well spotted. I made the fix :)

[dougiesquire]

LGTM - nice work!

[aekiss]

Great work both, LGTM!