"""Routines for calculating water fetch
LINCOM requires maps of fetch for calculating water roughness. Here we calculate this
using a fortran function, and include a utility function to create vector and raster
masks of water points based on a roughness raster map.
"""
import numpy as np
from windkit.metadata import update_var_attrs, _MAP_TYPE_ATTRS
from ..core import Rvea0330
from ..io.rastermap import rastermap_to_vectormap, rastermap_to_waspformat
from ..io.vectormap import vectormap_to_waspformat
from windkit.metadata import update_var_attrs, _MAP_TYPE_ATTRS
from .._errors import CoreError, _handle_fort_errors
FORT_FETCHMAP = Rvea0330.externals_mod.ffvec
[docs]
def calculate_fetchmap(vector_mask, raster_mask, wind):
"""Calculates a fetchmap for a given domain and wind
.. warning::
This function is experimental and its signature may change.
Parameters
----------
vector_mask : geopandas.GeoDataFrame or WaspVectorMap
vectormap of landmask (water (0) and land (1))
raster_mask : xarray.DataArray or WaspRasterMap
rastermap of landmask (water (0) and land (1))
wind : obj
LINCOM wind object
Returns
-------
xarray.DataArray
rastermap of the fetch
"""
fetch_map = rastermap_to_waspformat(raster_mask)
vector_mask = vectormap_to_waspformat(vector_mask)
# Get bounding box
xmin, ymin, xmax, ymax = fetch_map.bbox.bounds()
# Replace values with those from
fetch_map.values, ierr = FORT_FETCHMAP(
vector_mask.gx,
vector_mask.gy,
fetch_map.values,
fetch_map.nx,
fetch_map.ny,
xmin,
ymin,
xmax,
ymax,
wind.latitude,
wind.type,
wind.sector,
wind.wind_speed,
wind.wind_height,
wind.z0,
wind.west_east,
wind.south_north,
)
if ierr == 1:
_handle_fort_errors(95)
if ierr != 0:
_handle_fort_errors(ierr)
return fetch_map._to_rastermap()
[docs]
def roughness_to_landmask(lc_map, lctable):
"""Convert a roughness map to a landmask map
All classes with z0 = 0.0 will be converted to water (0) points, while all other
roughness lengths will be considered land.
.. warning::
This function is experimental and its signature may change.
Parameters
----------
lc_map: xarray.DataArray
Landcover rastermap
lctable: LandCoverTable
Map from LandCoverID to roughnesses
Returns
-------
raster_mask: pw.RasterMap
RasterMap of water (0) and land (1) mask
vector_mask: pw.VectorMap
VectorMap of water (0) and land (1) mask
"""
# TODO: Add vectormap support
ids = np.unique(lc_map)
# Create copy of data filled with ones
raster_mask = lc_map.copy()
raster_mask.values = np.ones_like(lc_map)
raster_mask.name = "landmask"
raster_mask = update_var_attrs(raster_mask, _MAP_TYPE_ATTRS)
# Replace ids that have roughness equal to 0
for k, v in lctable.items():
if k in ids and v["z0"] == 0:
raster_mask.values[lc_map.values == k] = 0
# Create vector mask from rastermap mask
vector_mask, _ = rastermap_to_vectormap(raster_mask)
return raster_mask, vector_mask