isofit.utils.skyview

Functions

skyview(input, output_directory[, resolution, ...])	Applies sky view factor calculation for a given projected DEM or DSM. Much of this code was borrowed from ARS Topo-Calc.
horizon_worker(angle, dem, spacing, aspect, tan_slope, ...)	Each worker gets an angle and is sent to this function to compute horizons, and save to a compressed/scaled zarr file.
init_horizon(output_directory, angles, shape)	Creates a zarr dataset for the horizon computations.
save_horizon(h, angle, output_directory)	utility function to house metadata and information for saving horizon angles
load_horizon(output_directory, angle)
create_shadow_mask(input, output_directory[, ...])	Computes horizons at a specific geometry to create a binary shadow mask because nearby terrain
update_h_for_local_topo(h, aspect, azimuth, tan_slope)
load_input(input, resolution[, obs_or_loc, method])
gradient_d8(dem, dx, dy[, aspect_rad])	Calculate the slope and aspect for provided dem,
calc_slope(dz_dx, dz_dy)	Calculate the slope given the finite differences
aspect(dz_dx, dz_dy)	Calculate the aspect from the finite difference.
aspect_to_ipw_radians(a)	IPW defines aspect differently than most GIS programs
horizon(azimuth, dem, spacing[, par, n_cores])	Calculate horizon angles for one direction. Horizon angles
hor1d(z[, fwd])
horval(z, delta, h)
hor2d(z, delta[, fwd, par, n_cores])
_hor2d_chunk(zbuf, delta, fwd, row_start, row_end)
adjust_spacing(spacing, skew_angle)	Adjust the grid spacing if a skew angle is present
skew(arr, angle[, fwd, fill_min])	Skew the origin of successive lines by a specified angle
skew_transpose(dem, spacing, angle)	Skew and transpose the dem for the given angle.
transpose_skew(dem, spacing, angle)	Transpose, skew then transpose a dem for the
cli(debug_args, **kwargs)

Module Contents

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.

Parameters:

• method (str, optional) – 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.

• n_angles (int, optional) – 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.

• keep_horizon_files (bool, optional) – 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.

• logging_level (str, optional) – Logging verbosity level (default is “INFO”); similar to apply_oe.

• log_file (str or None, optional) – File path to write logs; similar to apply_oe.

• n_cores (int, optional) – 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.

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.

init_horizon(output_directory, angles, shape)

Creates a zarr dataset for the horizon computations.

save_horizon(h, angle, output_directory)

utility function to house metadata and information for saving horizon angles

load_horizon(output_directory, angle)

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.

update_h_for_local_topo(h, aspect, azimuth, tan_slope)

load_input(input, resolution, obs_or_loc=None, method='slope')

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) )

Parameters:

• dem – array of elevation values

• dx – cell size along the x axis

• dy – cell size along the y axis

• aspect_rad – turn the aspect from degrees to IPW radians

Returns:

slope in radians aspect in degrees or IPW radians

calc_slope(dz_dx, dz_dy)

Calculate the slope given the finite differences

Parameters:

• dz_dx – finite difference in the x direction

• dz_dy – finite difference in the y direction

Returns:

slope numpy array

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/

Parameters:

• dz_dx – finite difference in the x direction

• dz_dy – finite difference in the y direction

Returns

aspect in degrees clockwise from North

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

Parameters:

a – aspect in degrees from due North

Returns

a: aspect in radians from due South

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.

Parameters:

• direction (azimuth {float} -- find horizon's along this)

• elevations (dem {np.array2d} -- numpy array of dem)

• spacing (spacing {float} -- grid)

Returns:

hcos {np.array} – cosines of angles to the horizon

hor1d(z, fwd=True)

horval(z, delta, h)

hor2d(z, delta, fwd=True, par=False, n_cores=1)

_hor2d_chunk(zbuf, delta, fwd, row_start, row_end)

adjust_spacing(spacing, skew_angle)

Adjust the grid spacing if a skew angle is present

Parameters:

• spacing (spacing {float} -- grid)

• [degrees] (skew_angle {float} -- angle to adjust the spacing for)

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

Parameters:

• arr – array to skew

• angle – angle between -45 and 45 to skew by

• fwd – add skew to image if True, unskew image if False

• fill_min – While IPW skew says it fills with zeros, the output image is filled with the minimum value

Returns:

skewed array

skew_transpose(dem, spacing, angle)

Skew and transpose the dem for the given angle. Also calculate the new spacing given the skew.

Parameters:

• elevations (dem {array} -- numpy array of dem)

• spacing (spacing {float} -- grid)

• angle (angle {float} -- skew)

Returns:

t – skew and transpose array spacing – new spacing adjusted for angle

transpose_skew(dem, spacing, angle)

Transpose, skew then transpose a dem for the given angle. Also calculate the new spacing

Parameters:

• elevations (dem {array} -- numpy array of dem)

• spacing (spacing {float} -- grid)

• angle (angle {float} -- skew)

Returns:

t – skew and transpose array spacing – new spacing adjusted for angle

cli(debug_args, **kwargs)