Source code for pywasp.lincom.grid

"""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( != 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, return out_da