"""Stores gridded topographic-data as an xarray object"""
import pyproj
import xarray as xr
import numpy as np
import windkit as wk
LINCOM_MAP_TYPE = {"elevation": 0, "roughness": 1, "fetch": 2}
[docs]
def open_raster(mapfilename, srs=None):
"""Reads grid from raster file
.. warning::
This function is experimental and its signature may change.
Parameters
----------
mapfilename: Path
Path to the rasterfile to be opened
srs : int, dict, str, or pyproj.CRS
Value to create pyproj.CRS object or an existing pyproj.CRS object.
(Default: use the CRS from the file.)
Returns
-------
xr.DataArray
DataArray containing the topographic data of interest
Notes
----
Open a rasterio file and rename dimensions to match PyWAsP Requirements.
Also ensure that the raster is ordered from lower-left not upper-left.
"""
grid = xr.open_rasterio(mapfilename).squeeze().drop_vars("band")
grid = grid.rename({"x": "west_east", "y": "south_north"})
# Flip y-axis if negative transform
if grid.attrs["transform"][4] < 0:
grid = grid.sel(south_north=slice(None, None, -1))
trans = list(grid.attrs["transform"])
trans[4] *= -1
trans[-1] = float(grid.south_north[0]) - 0.5 * trans[4]
grid.attrs["transform"] = tuple(trans)
# Add CRS to object using CF conventions
if srs is None:
if "crs" not in grid.attrs:
err = "File does not have a CRS, so you must provide an srs."
raise ValueError(err)
else:
srs = grid.attrs["crs"]
else:
# Check that file crs matches provided
if "crs" in grid.attrs:
if pyproj.CRS.from_string(grid.crs) != pyproj.CRS.from_user_input(srs):
err = "Provided srs, does not match file CRS."
raise ValueError(err)
del grid.attr["crs"]
if srs.startswith("+init="):
srs = srs[6:]
wk.spatial.add_crs(grid, srs)
return grid
[docs]
def grid_from_wasp_rastermap(rastermap):
"""Convert a WAsP RasterMap object to a Lincom Gridmap xarray object
.. warning::
This function is experimental and its signature may change.
Parameters
----------
rastermap: pw.RasterMap
Rastermap to be converted
Returns
-------
xarray.DataArray
DataArray containing the topographic map
"""
# Create coordinate arrays
x_points = np.arange(rastermap.nx) * rastermap.res + rastermap.minx
west_east = xr.DataArray(x_points, dims=("west_east"))
y_points = np.arange(rastermap.ny) * rastermap.res + rastermap.miny
south_north = xr.DataArray(y_points, dims=("south_north"))
# Create simple variables
out_da = xr.DataArray(
data=np.reshape(rastermap.values, (rastermap.ny, rastermap.nx)),
dims=("south_north", "west_east"),
# name='roughness',
# attrs={
# 'long_name': 'Dynamic Roughness',
# 'units': 'm',
# 'standard_name': 'surface_roughness_length'
# },
attrs={
"res": (rastermap.res, rastermap.res),
},
coords={"west_east": west_east, "south_north": south_north},
)
wk.spatial.add_crs(out_da, rastermap.crs)
return out_da