Visibility Analysis

Visibility analysis estimates which buildings or areas are visible from a given observer point (or set of points) within a specific distance. This is useful in assessing visual accessibility, urban form, and perceptual exposure in public space.

---

The module supports multiple modes of analysis:

Accurate Method

Computes visibility using fine-grained raster-based methods. More accurate for local areas, but slower.

objectnat.get_visibility_accurate(point_from, obstacles, view_distance, return_max_view_dist=False)

Function to get accurate visibility from a given point to buildings within a given distance.

Parameters:

Name	Type	Description	Default
point_from	Point | GeoDataFrame	The point or GeoDataFrame with 1 point from which the line of sight is drawn. If Point is provided it should be in the same crs as obstacles.	required
obstacles	GeoDataFrame	A GeoDataFrame containing the geometry of the obstacles.	required
view_distance	float	The distance of view from the point.	required
return_max_view_dist	bool	If True, the max view distance is returned with view polygon in tuple.	False

Returns:

Type	Description
Polygon | GeoDataFrame | tuple[Polygon | GeoDataFrame, float]	A polygon representing the area of visibility from the given point or polygon with max view distance. if point_from was a GeoDataFrame, return GeoDataFrame with one feature, else Polygon.

Notes

If a quick result is important, consider using the get_visibility() function instead. However, please note that get_visibility() may provide less accurate results.

Source code in src\objectnat\methods\visibility\visibility_analysis.py

def get_visibility_accurate(
    point_from: Point | gpd.GeoDataFrame, obstacles: gpd.GeoDataFrame, view_distance, return_max_view_dist=False
) -> Polygon | gpd.GeoDataFrame | tuple[Polygon | gpd.GeoDataFrame, float]:
    """
    Function to get accurate visibility from a given point to buildings within a given distance.

    Parameters:
        point_from (Point | gpd.GeoDataFrame):
            The point or GeoDataFrame with 1 point from which the line of sight is drawn.
            If Point is provided it should be in the same crs as obstacles.
        obstacles (gpd.GeoDataFrame):
            A GeoDataFrame containing the geometry of the obstacles.
        view_distance (float):
            The distance of view from the point.
        return_max_view_dist (bool):
            If True, the max view distance is returned with view polygon in tuple.

    Returns:
        (Polygon | gpd.GeoDataFrame | tuple[Polygon | gpd.GeoDataFrame, float]):
            A polygon representing the area of visibility from the given point or polygon with max view distance.
            if point_from was a GeoDataFrame, return GeoDataFrame with one feature, else Polygon.

    Notes:
        If a quick result is important, consider using the `get_visibility()` function instead.
        However, please note that `get_visibility()` may provide less accurate results.
    """

    def find_furthest_point(point_from, view_polygon):
        try:
            res = round(max(Point(coords).distance(point_from) for coords in view_polygon.exterior.coords), 1)
        except Exception as e:
            print(view_polygon)
            raise e
        return res

    local_crs = None
    original_crs = None
    return_gdf = False
    if isinstance(point_from, gpd.GeoDataFrame):
        original_crs = point_from.crs
        return_gdf = True
        if len(obstacles) > 0:
            local_crs = obstacles.estimate_utm_crs()
        else:
            local_crs = point_from.estimate_utm_crs()
        obstacles = obstacles.to_crs(local_crs)
        point_from = point_from.to_crs(local_crs)
        if len(point_from) > 1:
            logger.warning(
                f"This method processes only single point. The GeoDataFrame contains {len(point_from)} points - "
                "only the first geometry will be used for isochrone calculation. "
            )
        point_from = point_from.iloc[0].geometry
    else:
        obstacles = obstacles.copy()
    if obstacles.contains(point_from).any():
        return Polygon()
    obstacles.reset_index(inplace=True, drop=True)
    point_buffer = point_from.buffer(view_distance, resolution=32)
    allowed_geom_types = ["MultiPolygon", "Polygon", "LineString", "MultiLineString"]
    obstacles = obstacles[obstacles.geom_type.isin(allowed_geom_types)]
    s = obstacles.intersects(point_buffer)
    obstacles_in_buffer = obstacles.loc[s[s].index].geometry

    buildings_lines_in_buffer = gpd.GeoSeries(
        pd.Series(
            obstacles_in_buffer.apply(polygons_to_multilinestring).explode(index_parts=False).apply(explode_linestring)
        ).explode()
    )

    buildings_lines_in_buffer = buildings_lines_in_buffer.loc[buildings_lines_in_buffer.intersects(point_buffer)]

    buildings_in_buffer_points = gpd.GeoSeries(
        [Point(line.coords[0]) for line in buildings_lines_in_buffer.geometry]
        + [Point(line.coords[-1]) for line in buildings_lines_in_buffer.geometry]
    )

    max_dist = max(view_distance, buildings_in_buffer_points.distance(point_from).max())
    polygons = []
    buildings_lines_in_buffer = gpd.GeoDataFrame(geometry=buildings_lines_in_buffer, crs=obstacles.crs).reset_index()
    logger.debug("Calculation vis polygon")
    while not buildings_lines_in_buffer.empty:
        gdf_sindex = buildings_lines_in_buffer.sindex
        # TODO check if 2 walls are nearest and use the widest angle between points
        nearest_wall_sind = gdf_sindex.nearest(point_from, return_all=False, max_distance=max_dist)
        nearest_wall = buildings_lines_in_buffer.loc[nearest_wall_sind[1]].iloc[0]
        wall_points = [Point(coords) for coords in nearest_wall.geometry.coords]

        # Calculate angles and sort by angle
        points_with_angle = sorted(
            [(pt, math.atan2(pt.y - point_from.y, pt.x - point_from.x)) for pt in wall_points], key=lambda x: x[1]
        )
        delta_angle = 2 * math.pi + points_with_angle[0][1] - points_with_angle[-1][1]
        if round(delta_angle, 10) == round(math.pi, 10):
            wall_b_centroid = obstacles_in_buffer.loc[nearest_wall["index"]].centroid
            p1 = get_point_from_a_thorough_b(point_from, points_with_angle[0][0], max_dist)
            p2 = get_point_from_a_thorough_b(point_from, points_with_angle[1][0], max_dist)
            polygon = LineString([p1, p2])
            polygon = polygon.buffer(
                distance=max_dist * point_side_of_line(polygon, wall_b_centroid), single_sided=True
            )
        else:
            if delta_angle > math.pi:
                delta_angle = 2 * math.pi - delta_angle
            a = math.sqrt((max_dist**2) * (1 + (math.tan(delta_angle / 2) ** 2)))
            p1 = get_point_from_a_thorough_b(point_from, points_with_angle[0][0], a)
            p2 = get_point_from_a_thorough_b(point_from, points_with_angle[-1][0], a)
            polygon = Polygon([points_with_angle[0][0], p1, p2, points_with_angle[1][0]])

        polygons.append(polygon)
        buildings_lines_in_buffer.drop(nearest_wall_sind[1], inplace=True)

        if not polygon.is_valid or polygon.area < 1:
            buildings_lines_in_buffer.reset_index(drop=True, inplace=True)
            continue

        lines_to_kick = buildings_lines_in_buffer.within(polygon)
        buildings_lines_in_buffer = buildings_lines_in_buffer.loc[~lines_to_kick]
        buildings_lines_in_buffer.reset_index(drop=True, inplace=True)
    logger.debug("Done calculating!")
    res = point_buffer.difference(unary_union(polygons + obstacles_in_buffer.to_list()))

    if isinstance(res, MultiPolygon):
        res = list(res.geoms)
        for polygon in res:
            if polygon.intersects(point_from):
                res = polygon
                break

    if return_gdf:
        res = gpd.GeoDataFrame(geometry=[res], crs=local_crs).to_crs(original_crs)

    if return_max_view_dist:
        return res, find_furthest_point(point_from, res)
    return res

---

Fast Approximate Method

Optimized for large datasets or large areas. Uses geometry simplifications and vector-based visibility.

objectnat.get_visibility(point_from, obstacles, view_distance, resolution=32)

Function to get a quick estimate of visibility from a given point to buildings within a given distance.

Parameters:

Name	Type	Description	Default
point_from	Point | GeoDataFrame	The point or GeoDataFrame with 1 point from which the line of sight is drawn. If Point is provided it should be in the same crs as obstacles.	required
obstacles	GeoDataFrame	A GeoDataFrame containing the geometry of the buildings.	required
view_distance	float	The distance of view from the point.	required
resolution	int)	Buffer resolution for more accuracy (may give result slower)	32

Returns:

Type	Description
Polygon | GeoDataFrame	A polygon representing the area of visibility from the given point. if point_from was a GeoDataFrame, return GeoDataFrame with one feature, else Polygon.

Notes

This function provides a quicker but less accurate result compared to get_visibility_accurate(). If accuracy is important, consider using get_visibility_accurate() instead.

Source code in src\objectnat\methods\visibility\visibility_analysis.py

def get_visibility(
    point_from: Point | gpd.GeoDataFrame, obstacles: gpd.GeoDataFrame, view_distance: float, resolution: int = 32
) -> Polygon | gpd.GeoDataFrame:
    """
    Function to get a quick estimate of visibility from a given point to buildings within a given distance.

    Parameters:
        point_from (Point | gpd.GeoDataFrame):
            The point or GeoDataFrame with 1 point from which the line of sight is drawn.
            If Point is provided it should be in the same crs as obstacles.
        obstacles (gpd.GeoDataFrame):
            A GeoDataFrame containing the geometry of the buildings.
        view_distance (float):
            The distance of view from the point.
        resolution (int) :
            Buffer resolution for more accuracy (may give result slower)

    Returns:
        (Polygon | gpd.GeoDataFrame):
            A polygon representing the area of visibility from the given point.
            if point_from was a GeoDataFrame, return GeoDataFrame with one feature, else Polygon.

    Notes:
        This function provides a quicker but less accurate result compared to `get_visibility_accurate()`.
        If accuracy is important, consider using `get_visibility_accurate()` instead.
    """
    return_gdf = False
    if isinstance(point_from, gpd.GeoDataFrame):
        original_crs = point_from.crs
        return_gdf = True
        if len(obstacles) > 0:
            local_crs = obstacles.estimate_utm_crs()
        else:
            local_crs = point_from.estimate_utm_crs()
        obstacles = obstacles.to_crs(local_crs)
        point_from = point_from.to_crs(local_crs)
        if len(point_from) > 1:
            logger.warning(
                f"This method processes only single point. The GeoDataFrame contains {len(point_from)} points - "
                "only the first geometry will be used for isochrone calculation. "
            )
        point_from = point_from.iloc[0].geometry
    else:
        obstacles = obstacles.copy()
    point_buffer = point_from.buffer(view_distance, resolution=resolution)
    s = obstacles.intersects(point_buffer)
    buildings_in_buffer = obstacles.loc[s[s].index].reset_index(drop=True)
    buffer_exterior_ = list(point_buffer.exterior.coords)
    line_geometry = [LineString([point_from, ext]) for ext in buffer_exterior_]
    buffer_lines_gdf = gpd.GeoDataFrame(geometry=line_geometry)
    united_buildings = buildings_in_buffer.union_all()
    if united_buildings:
        splited_lines = buffer_lines_gdf["geometry"].apply(lambda x: x.difference(united_buildings))
    else:
        splited_lines = buffer_lines_gdf["geometry"]

    splited_lines_gdf = gpd.GeoDataFrame(geometry=splited_lines).explode(index_parts=True)
    splited_lines_list = []

    for _, v in splited_lines_gdf.groupby(level=0):
        splited_lines_list.append(v.iloc[0]["geometry"].coords[-1])
    circuit = Polygon(splited_lines_list)
    if united_buildings:
        circuit = circuit.difference(united_buildings)

    if return_gdf:
        circuit = gpd.GeoDataFrame(geometry=[circuit], crs=local_crs).to_crs(original_crs)
    return circuit

---

Catchment Visibility from Multiple Points

Performs visibility analysis for a dense grid of observer points.
Used to generate catchment areas of visibility (e.g., “where can this building be seen from?”).

objectnat.get_visibilities_from_points(points, obstacles, view_distance, sectors_n=None, max_workers=cpu_count())

Calculate visibility polygons from a set of points considering obstacles within a specified view distance.

Parameters:

Name	Type	Description	Default
points	GeoDataFrame	GeoDataFrame containing the points from which visibility is calculated.	required
obstacles	GeoDataFrame	GeoDataFrame containing the obstacles that block visibility.	required
view_distance	int	The maximum distance from each point within which visibility is calculated.	required
sectors_n	int	Number of sectors to divide the view into for more detailed visibility calculations. Defaults to None.	None
max_workers	int	Maximum workers in multiproccesing, multipocessing.cpu_count() by default.	cpu_count()

Returns:

Type	Description
list[Polygon]	A list of visibility polygons for each input point.

Notes

This function uses get_visibility_accurate() in multiprocessing way.

Source code in src\objectnat\methods\visibility\visibility_analysis.py

def get_visibilities_from_points(
    points: gpd.GeoDataFrame,
    obstacles: gpd.GeoDataFrame,
    view_distance: int,
    sectors_n=None,
    max_workers: int = cpu_count(),
) -> list[Polygon]:
    """
    Calculate visibility polygons from a set of points considering obstacles within a specified view distance.

    Parameters:
        points (gpd.GeoDataFrame):
            GeoDataFrame containing the points from which visibility is calculated.
        obstacles (gpd.GeoDataFrame):
            GeoDataFrame containing the obstacles that block visibility.
        view_distance (int):
            The maximum distance from each point within which visibility is calculated.
        sectors_n (int, optional):
            Number of sectors to divide the view into for more detailed visibility calculations. Defaults to None.
        max_workers (int, optional):
            Maximum workers in multiproccesing, multipocessing.cpu_count() by default.

    Returns:
        (list[Polygon]):
            A list of visibility polygons for each input point.

    Notes:
        This function uses `get_visibility_accurate()` in multiprocessing way.

    """
    if points.crs != obstacles.crs:
        raise ValueError(f"CRS mismatch, points crs:{points.crs} != obstacles crs:{obstacles.crs}")
    if points.crs.is_geographic:
        logger.warning("Points crs is geographic, it may produce invalid results")
    # remove points inside polygons
    joined = gpd.sjoin(points, obstacles, how="left", predicate="intersects")
    points = joined[joined.index_right.isnull()]

    # remove unused obstacles
    points_view = points.geometry.buffer(view_distance).union_all()
    s = obstacles.intersects(points_view)
    buildings_in_buffer = obstacles.loc[s[s].index].reset_index(drop=True)

    buildings_in_buffer.geometry = buildings_in_buffer.geometry.apply(
        lambda geom: MultiPolygon([geom]) if isinstance(geom, Polygon) else geom
    )
    args = [(point, buildings_in_buffer, view_distance, sectors_n) for point in points.geometry]
    all_visions = process_map(
        _multiprocess_get_vis,
        args,
        chunksize=5,
        desc="Calculating Visibility Catchment Area from each Point, it might take a while for a "
        "big amount of points",
        max_workers=max_workers,
    )

    # could return sectorized visions if sectors_n is set
    return all_visions

---

The image below shows an example of using visibility polygons to calculate "visibility pools" - areas in an urban environment that are most visible from different locations.

---