[Feature]: Explore land sea mask generation #576
Description
Is your feature request related to a problem?
I know that we had discussed this before and concluded to take an independent path to host land sea mask generation capability in a separate package (i.e., pcmdi_utils), I recently came across much simpler and more reliable method that possibly be worth considering to implement to xcdat.
The new approach uses regionmask, a package that is available via conda and pip.
This method has following advantages:
(1) It overcomes the complexity of the regrid2-based method that is originated from cdutil. Although the more precised land sea fraction conservation is considered in the regrid2-based method, in common practical case of global uses its influence is not very critical.
(2) It is not using global-land-mask, which is only available via pip, thus does not complicate xcdat installation.
Describe the solution you'd like
Function
[Proposed function updated -- boolean option added]
import regionmask
import xarray as xr
import xcdat as xc
def create_land_sea_mask(ds: xr.Dataset, boolean: bool=False) -> xr.DataArray:
"""
A function generates land sea mask (1: land, 0: sea) for given xarray Dataset,
assuming the given xarray dataset and has latitude and longitude coordinates.
Parameters
----------
ds : xr.Dataset
A Dataset object.
boolen : bool, optional
Set mask value to True (land) or False (sea), by default False
Returns
-------
xr.DataArray
A DataArray of land sea mask (1: land, 0: sea)
"""
# Create a land-sea mask using regionmask
land_mask = regionmask.defined_regions.natural_earth_v5_0_0.land_110
# Get the longitude and latitude from the xarray dataset
key_lon = xc.axis.get_dim_keys(ds, axis="X")
key_lat = xc.axis.get_dim_keys(ds, axis="Y")
lon = ds[key_lon]
lat = ds[key_lat]
# Mask the land-sea mask to match the dataset's coordinates
land_sea_mask = land_mask.mask(lon, lat)
if not boolean:
# Convert the land-sea mask to a boolean mask
land_sea_mask = xr.where(land_sea_mask, 0, 1)
return land_sea_maskExamples
[Examples updated -- boolean option and regional use case added]
target_grid = xc.create_uniform_grid(-90, 90, 1, 0, 359, 1)
mask = create_land_sea_mask(target_grid)
mask.plot()
target_grid = xc.create_uniform_grid(-90, 90, 0.5, -180, 179, 0.5)
mask = create_land_sea_mask(target_grid)
mask.plot()
mask = create_land_sea_mask(target_grid, boolean=True)
mask.plot()
target_grid = xc.create_uniform_grid(20, 45, 0.1, 110, 135, 0.1)
mask = create_land_sea_mask(target_grid, boolean=True)
mask.plot()
Describe alternatives you've considered
No response
Additional context
No response
Metadata
Metadata
Assignees
Type
Projects
Status