Smoothing of the Land-Sea Mask Dataset

The land-sea mask created from the ETOPO dataset comprises of 1s and 0s, where 1 denotes land and 0 denotes the ocean. However, it may often be useful to sort grid points not only by their land-surface type, but how distant they are from the coastline. GeoRegions.jl uses Gaussian Filtering from ImageFiltering.jl to perform a smoothing of the land-sea mask.

The degree of smoothing is expressed by the arguments σlon and σlat, which are integer inputs respectively.

You can get a smoothed land-sea mask via two means:

1. Calling the smoothing directly when retrieving the Land-Sea Dataset (recommended)
2. Smoothing a preexisting Land-Sea Dataset that has been loaded (not recommended)

The API for smoothing is listed below. However, it is worth nothing that when implementing the smoothing, you need to first call a larger GeoRegion around the GeoRegion of interest as a buffer.

Setup

using GeoRegions
using DelimitedFiles
using CairoMakie

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

Example comparison between Smooth and Unsmoothed Masks

geo = RectRegion("ACH","GLB","Aceh",[7,2,99,94],savegeo=false)
lsd_raw = getLandSea(geo,savelsd=false)
lsd_σ05 = getLandSea(geo,savelsd=false,smooth=true,σlon=5,σlat=5)
lsd_σ10 = getLandSea(geo,savelsd=false,smooth=true,σlon=10,σlat=10)

The Land-Sea Dataset (with Topography) has the following properties:
    Longitude Points    (lon) : Float32[94.00833, 94.025, 94.041664, 94.058334, 94.075, 94.09167, 94.10833, 94.125, 94.14167, 94.15833  …  98.84167, 98.85833, 98.875, 98.89167, 98.90833, 98.925, 98.941666, 98.958336, 98.975, 98.99167]
    Latitude Points     (lat) : Float32[2.0083334, 2.025, 2.0416667, 2.0583334, 2.075, 2.0916667, 2.1083333, 2.125, 2.1416667, 2.1583333  …  6.8416667, 6.858333, 6.875, 6.891667, 6.9083333, 6.925, 6.9416666, 6.9583335, 6.975, 6.991667]
    Region Size (nlon * nlat) : 300 lon points x 300 lat points

fig = Figure()

ax1 = Axis(
    fig[1,1],width=225,height=225,
    title="Raw Land-Sea Mask",xlabel="Longitude / º",ylabel="Latitude / º",
    limits=(93.9,99.1,1.9,7.1)
)
contourf!(
    ax1,lsd_raw.lon,lsd_raw.lat,lsd_raw.lsm,
    levels=range(0.05,0.95,length=19),extendlow=:auto,extendhigh=:auto
)
lines!(ax1,clon,clat,color=:black,linewidth=2)

ax2 = Axis(
    fig[1,2],width=225,height=225,
    title="Smoothed (σ=5)",xlabel="Longitude / º",
    limits=(93.9,99.1,1.9,7.1)
)
contourf!(
    ax2,lsd_σ05.lon,lsd_σ05.lat,lsd_σ05.lsm,
    levels=range(0.05,0.95,length=19),extendlow=:auto,extendhigh=:auto
)
lines!(ax2,clon,clat,color=:black,linewidth=2)

ax3 = Axis(
    fig[1,3],width=225,height=225,
    title="Smoothed (σ=10)",xlabel="Longitude / º",
    limits=(93.9,99.1,1.9,7.1)
)
contourf!(
    ax3,lsd_σ10.lon,lsd_σ10.lat,lsd_σ10.lsm,
    levels=range(0.05,0.95,length=19),extendlow=:auto,extendhigh=:auto
)
lines!(ax3,clon,clat,color=:black,linewidth=2)

resize_to_layout!(fig)
fig

Smoothing Directly from a loaded Land-Sea Mask

smooth!(lsd_raw,σlon=5,σlat=5)

f2 = Figure()

ax1 = Axis(
    f2[1,1],width=350,height=350,
    title="Smoothed from Global",xlabel="Longitude / º",ylabel="Latitude / º",
    limits=(93.9,99.1,1.9,7.1)
)
contourf!(
    ax1,lsd_σ05.lon,lsd_σ05.lat,lsd_σ05.lsm,
    levels=range(0.05,0.95,length=19),extendlow=:auto,extendhigh=:auto
)
lines!(ax1,clon,clat,color=:black,linewidth=2)

ax2 = Axis(
    f2[1,2],width=350,height=350,
    title="Directly Smoothed",xlabel="Longitude / º",
    limits=(93.9,99.1,1.9,7.1)
)
contourf!(
    ax2,lsd_raw.lon,lsd_raw.lat,lsd_raw.lsm,
    levels=range(0.05,0.95,length=19),extendlow=:auto,extendhigh=:auto
)
lines!(ax2,clon,clat,color=:black,linewidth=2)

resize_to_layout!(f2)
f2

You see here that smoothing the land-sea mask directly upon itself causes errors at the edges of the grid. This is because the smoothing applied is a circular smoothing, meaning that the boundaries are considered to be doubly-periodic, and thus bleed into each other. Thus, it is always best practice to call the smoothing directly when retrieving the Land-Sea dataset.

API

smooth!(
    lsd  :: LandSeaTopo;
    σlon :: Int = 0,
    σlat :: Int = 0,
    iterations :: Int = 100,
    topography :: Bool = false
) --> nothing

smooth!(
    lsd  :: LandSeaFlat;
    σlon :: Int = 0,
    σlat :: Int = 0,
    iterations :: Int = 100
) --> nothing

smooth!(
    lsm  :: Array{<:Real,2};
    σlon :: Int = 0,
    σlat :: Int = 0,
    iterations :: Int = 100
) --> nothing

smooth!(
    lsm  :: Array{<:Real,2},
    oro  :: Array{<:Real,2};
    σlon :: Int = 0,
    σlat :: Int = 0,
    iterations :: Int = 100
) --> nothing

smoothlsm(
    oro  :: Array{<:Real,2};
    σlon :: Int = 0,
    σlat :: Int = 0,
    iterations :: Int = 100
) --> Array{<:Real,2}