Loading and Saving Land-Sea Mask Datasets

The Land-Sea Dataset can be obtained using the function getLandSea(). There are two options:

1. Saving the LandSea dataset as a local NetCDF file
2. Extract the LandSea dataset directly from the OPeNDAP servers

Option (1) has a longer initial cost, because you have to download the global file, and then use the data from the global file to extract out the information for the GeoRegion of interest. However, once the data has been downloaded once, you don't need to download the data again and this greatly speeds up latency time.

Option (2) however, is much faster for smaller GeoRegions when loading the data directly from OPeNDAP is much easier than downloading the entire dataset. However, loading data directly from OPeNDAP comes with its own noticeable latency, which can add up if done iteratively.

You can toggle between the two options using the keyword argument savelsd - true sets the function to option 1, and false is option 2.

See end of the page for the API.

For more details on smoothing the land-sea mask, such that distance from the actual coastline is better resolved, please see this page

Setup

using GeoRegions
using DelimitedFiles
using CairoMakie

download("https://raw.githubusercontent.com/natgeo-wong/GeoPlottingData/main/coastline_resl.txt","coast.cst")
coast = readdlm("coast.cst",comments=true)
clon  = coast[:,1]
clat  = coast[:,2]
nothing

Retrieving LandSea Example over Aceh

geo  = RectRegion("ACH","GLB","Aceh",[6,3,98,95],savegeo=false)
slon,slat = coordGeoRegion(geo)
lsd = getLandSea(geo,savelsd=false)

The Land-Sea Dataset (with Topography) has the following properties:
    Longitude Points    (lon) : Float32[95.00833, 95.025, 95.041664, 95.058334, 95.075, 95.09167, 95.10833, 95.125, 95.14167, 95.15833  …  97.84167, 97.85833, 97.875, 97.89167, 97.90833, 97.925, 97.941666, 97.958336, 97.975, 97.99167]
    Latitude Points     (lat) : Float32[3.0083334, 3.025, 3.0416667, 3.0583334, 3.075, 3.0916667, 3.1083333, 3.125, 3.1416667, 3.1583333  …  5.8416667, 5.858333, 5.875, 5.891667, 5.9083333, 5.925, 5.9416666, 5.9583335, 5.975, 5.991667]
    Region Size (nlon * nlat) : 180 lon points x 180 lat points

And we plot it here, with coarser-resolution coastlines for comparison:

fig = Figure()
aspect = (maximum(slon)-minimum(slon))/(maximum(slat)-minimum(slat))

ax1 = Axis(
    fig[1,1],width=350,height=350/aspect,
    title="Orographic Height",xlabel="Longitude / º",ylabel="Latitude / º",
    limits=(minimum(slon)-0.1,maximum(slon)+0.1,minimum(slat)-0.1,maximum(slat)+0.1)
)
contourf!(
    ax1,lsd.lon,lsd.lat,lsd.z,
    levels=range(0,1500,length=16),extendlow=:auto,extendhigh=:auto
)
lines!(ax1,clon,clat,color=:black,linewidth=2)
lines!(ax1,slon,slat,linewidth=5)

ax2 = Axis(
    fig[1,2],width=350,height=350/aspect,
    title="Land-Sea Mask",xlabel="Longitude / º",
    limits=(minimum(slon)-0.1,maximum(slon)+0.1,minimum(slat)-0.1,maximum(slat)+0.1)
)
contourf!(
    ax2,lsd.lon,lsd.lat,lsd.lsm,
    levels=range(0.05,0.95,length=19),extendlow=:auto,extendhigh=:auto
)
lines!(ax2,slon,slat,linewidth=5)

resize_to_layout!(fig)
fig

API

getLandSea(
    geo  :: GeoRegion = GeoRegion("GLB");
    path :: AbstractString = homedir(),
    resolution :: Int = 60,
    bedrock   :: Bool = false,
    geoid     :: Bool = false,
    returnlsd :: Bool = false,
    savelsd   :: Bool = false,
    smooth    :: Bool = false,
    σlon :: Int = 0,
    σlat :: Int = 0,
    iterations :: Int = 100,
    FT = Float32
) -> LandSea