JuliaClimate/ClimateBase.jl

[FR] Land Mask

Opened this issue · 8 comments

Datseris commented

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.

Datseris commented

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

Balinus commented

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 commented

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 commented

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 commented

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 commented

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 commented

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.

  • 🎉1