Set! discontinuous across prime meridian

[milankl]

orography = NoOrography(spectral_grid)
model = PrimitiveWetModel(;spectral_grid, orography)
H, λ₀, φ₀, σ = 4000, 0, 50, 15         # height, lon, lat position, and width
set!(model, orography=(λ,φ) -> H*exp((-(λ-λ₀)^2 - (φ-φ₀)^2)/2σ^2), add=true)
heatmap(model.orography.orography)

has a discontinuity at 0˚E

@maximilian-gelbrecht you think we can fix that within the set! function or does one need to account for that in the (λ,φ) -> function definition?

[maximilian-gelbrecht]

Huh, I'll look into it. There has to be a possibility to do that neatly within set!

[maximilian-gelbrecht]

Okay, after thinking about this a bit more. This is too poorly defined mathematically to do within set!. You'd have to provide more additional information of how you want to have the continuity to be enforced and this just becomes messy.

I'd just add a "ghost point" to the function definition in the individual cases

using SpeedyWeather 
spectral_grid = SpectralGrid()
orography = NoOrography(spectral_grid)
model = PrimitiveWetModel(;spectral_grid, orography)
H, λ₀, φ₀, σ = 4000, 0, 50, 15         # height, lon, lat position, and width
set!(model, orography=(λ,φ) -> H*(exp((-(λ-λ₀)^2 - (φ-φ₀)^2)/2σ^2) + exp((-(λ-(λ₀+360))^2 - (φ-φ₀)^2)/2σ^2)), add=true)
plot(model.orography.orography)

`

[maximilian-gelbrecht]

It's not easy to generalise this to arbitrary functions, so I'd prefer to keep this out of set! I think.

One thing we could think to implement is a function that shifts a RingGrid by a certain longitude. I am not even sure if this is easy to do with all reduced grids, but it would be helpful in cases in which you - in some way - want the prime meridian to be somewhere else to solve a discontinuity.

[milankl]

One thing we could think to implement is a function that shifts a RingGrid by a certain longitude. I am not even sure if this is easy to do with all reduced grids, but it would be helpful in cases in which you - in some way - want the prime meridian to be somewhere else to solve a discontinuity.

All reduced grids we have have 4 faces around the poles, so you could shift longitudinally by 90, 180, or 270˚, e.g.

Ring1: (1, 2, 3, 4) -> (4, 1, 2, 3)
Ring2: (5, 6, 7, 8, 9, 10, 11, 12) -> (11, 12, 5, 6, 7, 8, 9, 10)
Ring3: (13, ..., 24 ) -> (21, ..., 24, 13, ..., 20)

[treigerm]

Just putting my comment here from SpeedyWeather/RainMaker.jl#15 for persistence:

I think the problem here stems from the fact that the function (λ,φ) -> H*(exp((-(λ-λ₀)^2 - (φ-φ₀)^2)/2σ^2) + exp((-(λ-(λ₀+360))^2 - (φ-φ₀)^2)/2σ^2)) is using a Euclidean distance function when the underlying geometry of the space is actually non-Euclidean. So instead calculating the distances between (lon, lat) pairs using geodesics should fix this.

Using GeographicLib.jl

using GeographicLib
earth_geodesic((lon1, lat1), (lon2, lat2)) = GeographicLib.inverse(lon1, lat1, lon2, lat2)[:dist]

as the distance works. But this calculates the distance in meters so the width of the mountain needs to be adjusted

H, λ₀, φ₀, σ = 4000, 0, 50, 1_500_000
orography = NoOrography(spectral_grid)
model = PrimitiveWetModel(;spectral_grid, orography)
set!(model, orography=(λ,φ) -> H*exp((- earth_geodesic((λ,φ), mountain_loc)^2)/2σ^2), add=true) heatmap(model.orography.orography)

gives

[maximilian-gelbrecht]

Yes, that's the more proper way of doing that.

So, to answer Milan's initial question here: The functions should take care of that. GeographicLib.jl is a good find for that.

I guess we can close this issue?

[maximilian-gelbrecht]

I think with #588 we can close this now.