isofit.utils.skyview ==================== .. py:module:: isofit.utils.skyview Functions --------- .. autoapisummary:: isofit.utils.skyview.skyview isofit.utils.skyview.horizon_worker isofit.utils.skyview.init_horizon isofit.utils.skyview.save_horizon isofit.utils.skyview.load_horizon isofit.utils.skyview.create_shadow_mask isofit.utils.skyview.update_h_for_local_topo isofit.utils.skyview.load_input isofit.utils.skyview.gradient_d8 isofit.utils.skyview.calc_slope isofit.utils.skyview.aspect isofit.utils.skyview.aspect_to_ipw_radians isofit.utils.skyview.horizon isofit.utils.skyview.hor1d isofit.utils.skyview.horval isofit.utils.skyview.hor2d isofit.utils.skyview._hor2d_chunk isofit.utils.skyview.adjust_spacing isofit.utils.skyview.skew isofit.utils.skyview.skew_transpose isofit.utils.skyview.transpose_skew isofit.utils.skyview.cli Module Contents --------------- .. py:function:: skyview(input: str, output_directory: str, resolution: float = np.nan, obs_or_loc: str = None, method: str = 'slope', n_angles: int = 64, logging_level: str = 'INFO', log_file: str = None, n_cores: int = 1, keep_horizon_files: bool = False, ray_address: str = None, ray_redis_password: str = None, ray_temp_dir: str = None, ray_ip_head: str = None) Applies sky view factor calculation for a given projected DEM or DSM. Much of this code was borrowed from ARS Topo-Calc. The key thing here was to create a python-only, rasterio-free port of this that could be used within ISOFIT. We also included improvements that are current in Jeff Dozier's horizon method in Matlab (https://github.com/DozierJeff/Topographic-Horizons). Following suggestions from Dozier (2021), multiprocessing is leveraged here w.r.t. to n_angles rotating the image. As default, sky view is computed with n angles = 64 which in most cases is of sufficient accuracy to resolve but more angles may be used. Optionally to this horizon based method, one can pass method="slope" to compute a faster estimate that may be sufficent for regions with lower relief. The slope based estimate is simply, svf = cos^2(slope/2). Yet another option is to pass an ISOFIT "OBS" or "LOC" file as input and using the obs_or_loc arg. OBS files have slope data and can be used for method='slope' only. LOC files have elevation data and can be used for method='slope'. One can also use the full horizon method on the LOC file although this is not recommended because the edges miss information (a warning will be triggered in this case).  :param input: Path to the projected ENVI File. If `obs_or_loc` is None or "loc" input is elevation. If is "obs", then it's a slope product. :type input: str :param output_directory: Directory path for temporary files and outputs during processing; similar to apply_oe. :type output_directory: str :param resolution: Spatial resolution of the projected DEM in coordinate units (GSD in meters). Required for elevation input data. :type resolution: float, optional :param obs_or_loc: Options here are 'obs', 'loc', or None. Default is None. If 'obs' is selected, it will pick the slope data from index 6 in OBS file. If 'loc' is selected it well select the elevation data from index 2. None will assume a single band elevation data is passed. :type obs_or_loc: str, optional :param method: Options are either "horizon" or "slope". Passing "horizon" runs the full computation and is recommended for very steep terrain. Passing "slope"" runs the simplifed calculation of svf=cos^2(slope/2) and can be useful for more mild slopes. :type method: str, optional :param n_angles: Number of angles used in horizon calculations (default is 64). As a reference, n=72 computes every 5deg, n=64 every 5.6deg, n=32 every 11.25deg, etc. :type n_angles: int, optional :param keep_horizon_files: Horizon angles are created in output_dir as a zarr data structure. False deletes files, and True keeps them. These angles are based from zenith. :type keep_horizon_files: bool, optional :param logging_level: Logging verbosity level (default is "INFO"); similar to apply_oe. :type logging_level: str, optional :param log_file: File path to write logs; similar to apply_oe. :type log_file: str or None, optional :param n_cores: Number of CPU cores to use for parallel processing (default is 1). Only used for method="horizon". Note: n_cores should ideally not be larger than n_angles. :type n_cores: int, optional .. py:function:: horizon_worker(angle, dem, spacing, aspect, tan_slope, output_directory, logging_level, log_file) Each worker gets an angle and is sent to this function to compute horizons, and save to a compressed/scaled zarr file. .. py:function:: init_horizon(output_directory, angles, shape) Creates a zarr dataset for the horizon computations. .. py:function:: save_horizon(h, angle, output_directory) utility function to house metadata and information for saving horizon angles .. py:function:: load_horizon(output_directory, angle) .. py:function:: create_shadow_mask(input: str, output_directory: str, resolution: float = np.nan, sza: float = np.nan, saa: float = np.nan, logging_level: str = 'INFO', log_file: str = None, n_cores: int = 1, ray_address: str = None, ray_redis_password: str = None, ray_temp_dir: str = None, ray_ip_head: str = None) Computes horizons at a specific geometry to create a binary shadow mask because nearby terrain can cast shadows onto pixels as a function of the solar geometry that aren't always captured in cos-i. In this case, the input angle is the solar azimuth, and the solar zenith is compared to h at each pixel. As of right now, this assumes solar azimuth is constant value (0-360deg). However, solar zenith (0-90deg) can vary by pixel. .. py:function:: update_h_for_local_topo(h, aspect, azimuth, tan_slope) .. py:function:: load_input(input, resolution, obs_or_loc=None, method='slope') .. py:function:: gradient_d8(dem, dx, dy, aspect_rad=False) Calculate the slope and aspect for provided dem, using a 3x3 cell around the center Given a center cell e and it's neighbors: | a | b | c | | d | e | f | | g | h | i | The rate of change in the x direction is [dz/dx] = ((c + 2f + i) - (a + 2d + g) / (8 * dx) The rate of change in the y direction is [dz/dy] = ((g + 2h + i) - (a + 2b + c)) / (8 * dy) The slope is calculated slope_radians = arctan ( sqrt ([dz/dx]^2 + [dz/dy]^2) ) :param dem: array of elevation values :param dx: cell size along the x axis :param dy: cell size along the y axis :param aspect_rad: turn the aspect from degrees to IPW radians :returns: slope in radians aspect in degrees or IPW radians .. py:function:: calc_slope(dz_dx, dz_dy) Calculate the slope given the finite differences :param dz_dx: finite difference in the x direction :param dz_dy: finite difference in the y direction :returns: slope numpy array .. py:function:: aspect(dz_dx, dz_dy) Calculate the aspect from the finite difference. Aspect is degrees clockwise from North (0/360 degrees) See below for a referance to how ArcGIS calculates slope http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#/How_Aspect_works/00q900000023000000/ :param dz_dx: finite difference in the x direction :param dz_dy: finite difference in the y direction Returns aspect in degrees clockwise from North .. py:function:: aspect_to_ipw_radians(a) IPW defines aspect differently than most GIS programs so convert an aspect in degrees from due North (0/360) to the IPW definition. Aspect is radians from south (aspect 0 is toward the south) with range from -pi to pi, with negative values to the west and positive values to the east :param a: aspect in degrees from due North Returns a: aspect in radians from due South .. py:function:: horizon(azimuth, dem, spacing, par=False, n_cores=1) Calculate horizon angles for one direction. Horizon angles are based on Dozier and Frew 1990 and are adapted from the IPW C code. The coordinate system for the azimuth is 0 degrees is South, with positive angles through East and negative values through West. Azimuth values must be on the -180 -> 0 -> 180 range. :param azimuth {float} -- find horizon's along this direction: :param dem {np.array2d} -- numpy array of dem elevations: :param spacing {float} -- grid spacing: :returns: hcos {np.array} -- cosines of angles to the horizon .. py:function:: hor1d(z, fwd=True) .. py:function:: horval(z, delta, h) .. py:function:: hor2d(z, delta, fwd=True, par=False, n_cores=1) .. py:function:: _hor2d_chunk(zbuf, delta, fwd, row_start, row_end) .. py:function:: adjust_spacing(spacing, skew_angle) Adjust the grid spacing if a skew angle is present :param spacing {float} -- grid spacing: :param skew_angle {float} -- angle to adjust the spacing for [degrees]: .. py:function:: skew(arr, angle, fwd=True, fill_min=True) Skew the origin of successive lines by a specified angle A skew with angle of 30 degrees causes the following transformation: +-----------+ +---------------+ | | |000/ /| | input | |00/ output /0| | image | |0/ image /00| | | |/ /000| +-----------+ +---------------+ Calling skew with fwd=False will return the output image back to the input image. Skew angle must be between -45 and 45 degrees :param arr: array to skew :param angle: angle between -45 and 45 to skew by :param fwd: add skew to image if True, unskew image if False :param fill_min: While IPW skew says it fills with zeros, the output image is filled with the minimum value :returns: skewed array .. py:function:: skew_transpose(dem, spacing, angle) Skew and transpose the dem for the given angle. Also calculate the new spacing given the skew. :param dem {array} -- numpy array of dem elevations: :param spacing {float} -- grid spacing: :param angle {float} -- skew angle: :returns: t -- skew and transpose array spacing -- new spacing adjusted for angle .. py:function:: transpose_skew(dem, spacing, angle) Transpose, skew then transpose a dem for the given angle. Also calculate the new spacing :param dem {array} -- numpy array of dem elevations: :param spacing {float} -- grid spacing: :param angle {float} -- skew angle: :returns: t -- skew and transpose array spacing -- new spacing adjusted for angle .. py:function:: cli(debug_args, **kwargs)