"""
Contains class and class methods to manipulate elevation and roughness vector files
VectorMaps are used to represent elevation contours and roughness change lines.
Pywasp can read elevation and roughness data from either classic WAsP mapfiles
or the newer gml file formats, but many other file formats are also supported
through GDAL's OGR tool, however this requires you to set the field ids
appropriately when reading the files.
VectorMaps are used by pywaps Topography module to calculate speedups, turning,
and the mesoroughness, for either generalization or downscaling.
Mapfiles are used to represent elevation contours and roughness change lines.
They can then be used by the Topography module to calculate speedups, turning,
and the mesoroughness, for either generalization or downscaling.
"""
import warnings
import numpy as np
from .landcover import LandCoverTable
from ..core import Rvea0287, Rvea0290 #
from windkit.spatial import BBox
import pyproj
from shapely.geometry import LineString, Point, Polygon
from shapely import linestrings
from windkit._vectormap_helpers import (
_is_vectormap,
_get_map_type,
_z0_to_landcover,
explode_gdf,
)
from windkit._vectormap_helpers import (
_LR_COLS,
VECTORMAP_ELEV_COL,
_MAP_TYPE_CODES,
VECTORMAP_ROUL_COL,
VECTORMAP_ROUR_COL,
VECTORMAP_IDL_COL,
VECTORMAP_IDR_COL,
VECTORMAP_LMASKL_COL,
VECTORMAP_LMASKR_COL,
)
import logging
from inspect import cleandoc
from collections import defaultdict
import geopandas as gpd
[docs]
def vectormap_to_rastermap(
gdf, res, xmin=None, ymin=None, xmax=None, ymax=None, lctable=None
):
"""
Converts a geopandas.GeoDataFrame vectormap object to a xarray.DataArray
rastermap object.
.. warning::
This function is experimental and its signature may change.
Parameters
----------
gdf : geopandas.GeoDataFrame
vectormap to be coverted
res : float
Resolution of x and y dimensions (dx, dy)
xmin, ymin, xmax, ymax : float, optional
Bounding box values of the output raster map. By default (None), the
entire vectormap bounding box is used.
lctable : LandCoverTable
Landcover table specifying the landcover classes and their mappings to
roughness and displacement height. Required for 'landcover' maps.
Returns
-------
xarray.DataArray
Rastermap representation of the map
"""
wvm = vectormap_to_waspformat(gdf, lctable)
wrm = wvm.to_wasp_rastermap(res, xmin, ymin, xmax, ymax)
return wrm._to_rastermap()
class WaspVectorMap(object):
"""
Vector representation of the map.
Notes
-----
This is the only type of mapfile that the Topography module works with for
the moment.
Parameters
----------
header : str
Single line description of the mapfile, will typically be the proj.4
string that describes the projection of the map.
gx : array_like
The first three entries of this array are the left roughness length,
right roughness length, and elevation, followed by the x-cooridinate
values for each point. (1D)
gy : array_like
The first three entries of this array are the number of points for a
specific line, the line type, and a null value represented as 0.0. (1D)
n_lines : int
Number of lines in the map
n_pts : int
Total number of points in the map
epsg : int
EPSG number for identifying SRS
map_type : int
ID representing the type of map from MapType
bbox : tuple
Bounding box holding the global extent of the object. This has the
following values (minx, miny, maxx, maxy).
lctable : LandCoverTable
Landcover table that maps landcover to roughness
"""
_map_type_codes = {
"elevation": 0,
"roughness": 1,
"landcover": 6,
"landmask": 15,
}
def __init__(
self, header, gx, gy, n_lines, n_pts, crs, map_type, bbox, lctable=None
):
# Always store gx and gy as a list
if isinstance(gx, np.ndarray):
gx = gx.tolist()
if isinstance(gy, np.ndarray):
gy = gy.tolist()
if not isinstance(bbox, BBox):
raise ValueError("bbox value should be a windkit.BBox object")
self.header = header
self.gx = gx
self.gy = gy
self.n_lines = n_lines
self.n_pts = n_pts
self.crs = crs
if crs is not None:
self.crs = pyproj.CRS.from_user_input(crs)
self.map_type = map_type
self.bbox = bbox
self.lctable = lctable
def __repr__(self): # pragma: no cover
repr_list = [
repr(self.header),
repr(self.gx),
repr(self.gy),
repr(self.n_lines),
repr(self.n_pts),
repr(self.crs.to_string()),
repr(self.map_type),
repr(self.bbox),
]
return "WaspVectorMap(" + ", ".join(repr_list) + ")"
def __str__(self): # pragma: no cover
out = """
{0.header}
Contains {0.n_lines} contour lines
It has the extent of {0.bbox}
""".format(
self
)
return out
@classmethod
def from_geodf(cls, gdf, lctable=None):
"""
Initialize a WaspVectorMap from a geodataframe
.. warning::
This function is experimental and its signature may change.
Parameters
----------
gdf: geopandas.GeoDataFrame
Vectormap as a geodataframe
map_type: str
Either "elevation" or "roughness" to split combo map
"""
args, kwargs = _vectormap_to_waspformat(gdf, lctable)
return cls(*args, **kwargs)
def _to_geodf(self):
"""
Convert WaspVectorMap to GeoDataFrame object.
"""
return _waspformat_to_vectormap(
self.gx, self.gy, self.n_lines, self.crs, self.lctable
)
def to_wasp_rastermap(self, res, xmin=None, ymin=None, xmax=None, ymax=None):
"""
Converts a WaspVectorMap object to a WaspRasterMap object.
.. warning::
This function is experimental and its signature may change.
Parameters
----------
res : float
Resolution of x and y dimensions (dx, dy)
xmin, ymin, xmax, ymax : float, optional
Bounding box values of the output raster map. By default (None), the
entire WaspVectorMap bounding box is used.
Returns
-------
WaspRasterMap
RasterMap class instance
"""
from .rastermap import WaspRasterMap # Here to prevent circular imports
# if no bounding box is specified, just plot the whole map based on bbox from vectormap
if any([x is None for x in [xmin, ymin, xmax, ymax]]):
if not all([x is None for x in [xmin, ymin, xmax, ymax]]):
raise ValueError("Need to set all bounds if you are setting any.")
xmin, ymin, xmax, ymax = self.bbox.bounds()
nx = int((xmax - xmin) / res)
ny = int((ymax - ymin) / res)
if self.map_type == "elevation":
convert = Rvea0290.elegridder_mod.elevationmaptogrid
elif self.map_type in ["roughness", "landcover", "landmask"]:
convert = Rvea0290.rghgridder_mod.roughnessmaptogrid
else:
err = f"Don't know how to convert WaspVectorMap of map_type {self.map_type} to WaspRasterMap"
raise ValueError(err)
values, _ = convert(self.gx, self.gy, self.n_lines, xmin, ymin, nx, ny, res)
return WaspRasterMap(
values,
xmin,
ymin,
nx,
ny,
res,
self.crs,
self.map_type,
lctable=self.lctable,
)
def _to_landcovermap(self, lctable=None):
"""
Generate landcover map and table from WaspVectormap
Parameters
----------
lctable : dict
Dictionary representing lookup table
Returns
-------
lcmap : :any:`WaspVectorMap`
WaspVectorMap class instance with a defined LandCoverTable, by default None
Notes
-----
This function creates the landcover classes from a vector map by
looping over all vector lines and checking that the roughness
length from each line is represented in the resulting LandCoverTable.
Each of the roughness lines is replaced with a corresponding
landcover id, which is then tied to a row in the landcover table.
"""
warnings.warn(
"`_to_landcovermap` is deprecated. Use `_z0_to_landcover` from windkit instead!",
FutureWarning,
)
if lctable is not None and self.map_type == "landcover":
warnings.warn(
"""This is already a landcover
map, so will not convert anything"""
)
if lctable is not None and self.map_type == "roughness":
warnings.warn(
"""There is a LandCoverTable specified but
the map_type indicates it is a roughness map. The lctable
will be overwritten by the values that are in the
mapfile"""
)
if self.map_type != "roughness" and self.map_type != "landcover":
raise ValueError(
"""This is not a roughness map because
it has maptype {0}""".format(
self.map_type
)
)
gx = np.array(self.gx, dtype="f", order="F")
gy = np.array(self.gy, dtype="f", order="F")
if lctable is None and self.map_type == "roughness":
nlucs, table, ierr = Rvea0287.rmap_to_luc(gx, gy, self.n_lines)
if ierr == 1:
warnings.warn(
"""The number of maximum landcover classes in
the map was reached. Therefore it is
possible that not all roughness values
in the map have a lookup roughness in the
generated landcover table. Please make sure
that your map has less than 100
unique roughness lengths"""
)
diclctable = {}
for i in range(nlucs):
diclctable[table[i, 0]] = {
"z0": table[i, 1],
"d": table[i, 2],
"desc": "",
}
lctable = LandCoverTable(diclctable)
# return vector map with landcover id's instead of roughness length
lumap = WaspVectorMap(
self.header,
gx,
gy,
self.n_lines,
self.n_pts,
self.crs,
"landcover",
self.bbox,
lctable=lctable,
)
return lumap
def _vectormap_to_waspformat(gdf, lctable=None):
"""
Create arguments and keyword-agument
to create a WaspVectorMap from a geodataframe
Parameters
----------
gdf : geopandas.GeoDataFrame
Vectormap as a geodataframe
elev_col : str
name of elevation column in geodataframe
z0_cols : bool
names of z0_left and z0_right in geodataframe
Returns
-------
args : tuple
Arguments needed for WaspVectorMap __init__
kwargs : dict
Key-word aguments for WaspVectorMap __init__
"""
if not _is_vectormap(gdf):
raise ValueError(
cleandoc(
"""This is not a GeoDataFrame. Perhaps you are passing in a combination of a GeoDataFrame and a landcover table?"""
)
)
# Wasp only works with single part geometries
gdf = explode_gdf(gdf)
# Get the map_type based con column names, and throw errors / warnings based on type
map_type = _get_map_type(gdf)
if map_type == "roughness":
if lctable is not None:
logging.warning(
cleandoc(
"""You specified a landcover table with a roughness map, but the roughnesses specified in the roughness map will be used instead."""
)
)
elif map_type == "landcover":
if lctable is None:
raise ValueError(
"You need to specify a landcover table with a landcover map"
)
# If we have non-elevation maps we have left-right columns, get their names
if map_type in ["roughness", "landcover", "landmask"]:
lr_cols = _LR_COLS[map_type]
else:
lr_cols = 2 * (VECTORMAP_ELEV_COL,)
coords = gdf.get_coordinates()
gx, gy, _, n_lines, ierr = Rvea0287.waspvectormap.from_gpd(
coords.index,
coords.x.values,
coords.y.values,
gdf[lr_cols[0]].values,
gdf[lr_cols[1]].values,
_MAP_TYPE_CODES[map_type],
)
n_pts = coords.index.size
if ierr > 0:
raise ValueError("Unknown map type!")
crs = gdf.crs
bbox = BBox.from_cornerpts(*gdf.total_bounds, crs)
args = ("PyWAsP map", gx, gy, n_lines, n_pts, crs, map_type, bbox)
kwargs = {"lctable": lctable if map_type == "landcover" else None}
return args, kwargs
def _waspformat_to_vectormap(gx, gy, n_lines, crs=None, lctable=None):
"""
Create geopandas.GeoDataFrame from WaspVectorMap arrays
gx, gy, and n_lines (the number of features).
Parameters
----------
gx, gy : array_like
WaspVectorMap arrays
n_lines : int
Number of features
crs : str, int, tuple, dict, pyproj.CRS
Coordinate Reference System.
Returns
-------
geopandas.GeoDataFrame
Vector map in geopandas format.
"""
(
index,
x,
y,
z0_left,
z0_right,
elev,
map_type,
ierror,
) = Rvea0287.waspvectormap.to_gpd(gx, gy, len(gx) - n_lines * 3, n_lines)
if ierror == 1:
raise ValueError("""No WAsP lines found!""")
geom = linestrings(x, y, indices=index)
if all(map_type == 0):
data = {VECTORMAP_ELEV_COL: elev}
elif all(map_type == 1):
data = {VECTORMAP_ROUL_COL: z0_left, VECTORMAP_ROUR_COL: z0_right}
elif all(map_type == 6):
if lctable is None:
raise ValueError(
"You need to specify a landcover table with a landcover map"
)
data = {VECTORMAP_IDL_COL: z0_left, VECTORMAP_IDR_COL: z0_right}
elif all(map_type == 15):
data = {VECTORMAP_LMASKL_COL: z0_left, VECTORMAP_LMASKR_COL: z0_right}
else:
raise ValueError("More than 1 map type found!")
gdf = gpd.GeoDataFrame(data, geometry=geom, crs=crs)
if all(map_type == 0):
return gdf
elif all(map_type == 1):
return _z0_to_landcover(gdf)
return (gdf, lctable)
def waspformat_to_vectormap(wvm, lctable=None):
"""
Convert WaspVectorMap to geopandas.GeoDataFrame
.. warning::
This function is experimental and its signature may change.
Parameters
----------
wvm : WaspVectorMap
VectorMap to be converted to GeoDataFrame
Returns
-------
geopandas.GeoDataFrame
Resulting geopandas.GeoDataFrame
"""
return _waspformat_to_vectormap(
wvm.gx, wvm.gy, wvm.n_lines, crs=wvm.crs, lctable=lctable
)