from pathlib import Path
import geopandas as gpd
import numpy as np
from shapely import box
from shapely import contains
from shapely import disjoint
from shapely import points
[docs]
def deposition_coeff_from_land_occupation_data(path_data, file_name_data, land_occupation_to_deposition_coeff, msh):
r"""
Derive the map of the deposition coefficient :math:`\tau_{loss}(x,y)`
given the land occupation data
and the value of deposition coefficient with respect to the land occupation.
Read the land occupation data,
search the cells of the mesh contained in each polygon of the land occupation data
and return the associated deposition coefficient map.
Parameters
----------
path_data: str
Path to the folder that contains the land occupation data.
file_name_data: str
Name of the file that contains the land occupation data.
land_occupation_to_deposition_coeff: dict
Dictionary containing the deposition coefficient with respect to the land occupation.
msh: ~pheromone_dispersion.geom.MeshRect2D
The geometry of the domain.
Returns
-------
~numpy.ndarray
The value of the deposition coefficient :math:`\tau_{loss}(x,y)` at each cells of the mesh.
Raises
------
ValueError
if one of the land occupation has no deposition coefficient value,
i.e. one of the land occupation keyword of the land occupation data
is not a keyword of the dictionnary containing
the deposition coefficient with respect to the land occupation as well
Notes
-----
- Ensure that points are indeed in one of the polygons.
"""
# To do:
# Take into account the case None.
# Add exceptions to ensure points are in one of the polygons.
# Add exceptions to check the format of the dataframe.
# Add exceptions to verify that the land occupation data and the dictionary coincide.
df = gpd.read_file(Path(path_data) / file_name_data)
# removing the None, the associated points will then be setted to 0
# better way to do it ?
for i in df.index:
if df['CLAS_SCOT'][i] is None:
df = df.drop(i)
xm, ym = np.meshgrid(msh.x, msh.y)
geoms = points(xm, ym)
msh_box = box(np.min(msh.x), np.min(msh.y), np.max(msh.x), np.max(msh.y))
deposition_coeff = np.zeros((msh.y.size, msh.x.size), dtype=float)
for polygon, land_occupation in zip(df['geometry'], df['CLAS_SCOT']):
if not disjoint(msh_box, polygon):
if land_occupation not in land_occupation_to_deposition_coeff.keys():
raise ValueError(
"The land occupation classification "
+ str(land_occupation)
+ "is not found in the dictionary containing the deposition coefficient with respect to the land occupation"
)
deposition_coeff[contains(polygon, geoms)] = land_occupation_to_deposition_coeff[str(land_occupation)]
return deposition_coeff