Skip to content

Skyview

source module isofit.utils.skyview

Functions

  • skyview — 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.

  • horizon_worker — Each worker gets an angle and is sent to this function to compute horizons, and save to a compressed/scaled netcdf file.

  • save_horizon — utility function to house metadata and information for saving horizon angles

  • load_horizon — load horizon netcdf in for skyview calc using scale factors defined in save.

  • create_shadow_mask — 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.

  • update_h_for_local_topo

  • load_input

  • gradient_d8 — Calculate the slope and aspect for provided dem, using a 3x3 cell around the center

  • calc_slope — Calculate the slope given the finite differences

  • aspect — Calculate the aspect from the finite difference. Aspect is degrees clockwise from North (0/360 degrees)

  • aspect_to_ipw_radians — IPW defines aspect differently than most GIS programs so convert an aspect in degrees from due North (0/360) to the IPW definition.

  • horizon — Calculate horizon angles for one direction. Horizon angles are based on Dozier and Frew 1990 and are adapted from the IPW C code.

  • hor1d

  • horval

  • hor2d

  • adjust_spacing — Adjust the grid spacing if a skew angle is present

  • skew — Skew the origin of successive lines by a specified angle A skew with angle of 30 degrees causes the following transformation:

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

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

  • cli

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

 Parameters

input : str 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. output_directory : str Directory path for temporary files and outputs during processing; similar to apply_oe. resolution : float, optional Spatial resolution of the projected DEM in coordinate units (GSD in meters). Required for elevation input data. obs_or_loc : str, optional 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. 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 netcdfs. 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.

Raises

  • TypeError

  • ValueError

source 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 netcdf file.

source save_horizon(h, angle, output_directory)

utility function to house metadata and information for saving horizon angles

source load_horizon(file_path)

load horizon netcdf in for skyview calc using scale factors defined in save.

source 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.

Raises

  • ValueError

source update_h_for_local_topo(h, aspect, azimuth, tan_slope)

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

Raises

  • FileNotFoundError

  • ValueError

source 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

source 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

source 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

source 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

source 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

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

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

  • spacing {float} -- grid spacing

Returns

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

Raises

  • ValueError

source hor1d(z, fwd=True)

source horval(z, delta, h)

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

Raises

  • ValueError

source adjust_spacing(spacing, skew_angle)

Adjust the grid spacing if a skew angle is present

Parameters

  • spacing {float} -- grid spacing

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

Raises

  • ValueError

source 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

Raises

  • ValueError

source skew_transpose(dem, spacing, angle)

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

Parameters

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

  • spacing {float} -- grid spacing

  • angle {float} -- skew angle

Returns

  • t -- skew and transpose array spacing -- new spacing adjusted for angle

source transpose_skew(dem, spacing, angle)

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

Parameters

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

  • spacing {float} -- grid spacing

  • angle {float} -- skew angle

Returns

  • t -- skew and transpose array spacing -- new spacing adjusted for angle

source cli(debug_args, **kwargs)