[FR] Land Mask

[Datseris]

With current functionality it is straight-forward to make a "land mask". This means to average given array R over only land or over only ocean. To do this we first create a filed W that is 0 over ocean and 1 over land.

W = [island(lon, lat) ? 1.0 : 0.0 for lon in lons for lat in lats]
Rland = spacemean(R, W)

What we actually need is to create a function island that returns 1 if the given longitude and latitude is over land, and 0 otherwise. How to do it?

The alternative would be to save a high resolution lon-lat grid with booleans to disk, and load this every time and check. But I think this file might be too large to be attached to the source code of ClimateBase.jl.

---

So, depending on how easy it is to make the island function, this issue is either easy or hard.

[Balinus]

Some dataset about land classes : http://www.fao.org/land-water/land/land-governance/land-resources-planning-toolbox/category/details/en/c/1036354/

[Datseris]

My concern is, should we bundle a really large .nc dataset with this package? Then upon install users would have to download it...

[Balinus]

Could it be downloaded on a need basis? For example, cartopy download land-mask when called the 1st time. We'll probably just need a message when that happen. How big is the file? I have not seen it.

[felixcremer]

We might use the https://github.com/JuliaGeo/GeoDatasets.jl package to get a land sea mask. But I am not sure, whether this is fine grained enough for this purpose.

[Datseris]

https://juliageo.org/GeoDatasets.jl/dev/#GeoDatasets.landseamask , given in arc minutes

Yeah, it is definitely fine grained enough at least for my scope. @Balinus what do you think?

[Balinus]

Seems interesting! Never heard of GeoDatasets. arc minute seems high resolution enough for the purpose of climate modelling (at least, until convection-permitting models becomes more available).

As a side note, the Julia ecosystem is really flourishing in the last couple of months. There is so much new packages, I can't keep up.

[Datseris]

As a side note, the Julia ecosystem is really flourishing in the last couple of months. There is so much new packages, I can't keep up.

Absolutely agreed, and I share the sentiment, especially given that I have to spend most of my open-source related time on JuliaDynamics :( E.g. bridging with the work of Meshes / MeshGrids and GeoStats would be amazing, but it feels such a daunting task for me at the moment!!!

[Datseris]

So, we have found an easy way to download "large" files on demand when a function is used, without bundling them with Pkg.add, using artifacts. We did it in e.g., Agents.jl:

https://github.com/JuliaDynamics/Agents.jl/blob/1f6eb1d43c5776ced6791af7cbb54a4d72a590cb/src/spaces/openstreetmap.jl#L116-L137

"""
    OSM.test_map()

Download a small test map of [Göttingen](https://nominatim.openstreetmap.org/ui/details.html?osmtype=R&osmid=191361&class=boundary)
as an artifact. Return a path to the downloaded file.
"""
function test_map()
    artifact_toml = joinpath(@__DIR__, "../../Artifacts.toml")
    map_hash = artifact_hash("osm_map_gottingen", artifact_toml)
    if isnothing(map_hash) || !artifact_exists(map_hash)
        map_hash = create_artifact() do artifact_dir
            Downloads.download(
                "https://raw.githubusercontent.com/JuliaDynamics/JuliaDynamics/master/artifacts/agents/osm_map_gottingen.json",
                joinpath(artifact_dir, "osm_map_gottingen.json")
            )
        end


        bind_artifact!(artifact_toml, "osm_map_gottingen", map_hash; force = true)
    end


    return joinpath(artifact_path(map_hash), "osm_map_gottingen.json")
end

This seems to be the natural way to implement a land mask. We could also have several land masks stored, with different resultions.