შიგთავსზე გადასვლა

👁️ ხილვადობა (Viewshed)

ხედვის არე (viewshed) არის რუკის ყველა ის ადგილი, რომელიც მოცემული წერტილიდან ჩანს. ამას იყენებენ სმარტფონის ანძების განლაგებაში, ხანძარსაწინააღმდეგო კოშკურების დაგეგმვაში, ქარის ტურბინების ვიზუალური გავლენის შეფასებაში და ტურისტული ხედების ანალიზში. ყველაფერი ერთ კითხვაზე დგას: ხედავს თუ არა წერტილი A წერტილ B-ს, თუ რელიეფი ხელს უშლის?

💡 იდეა: ხედვის ხაზი (line of sight)

A-დან B-მდე გავავლოთ სწორი ხაზი და ვნახოთ, რელიეფი რომელიმე წერტილში ამ ხაზს ზემოთ ხომ არ ადის. ამისთვის რელიეფის ყოველ უჯრაზე ვითვლით დახრილობას დამკვირვებლის თვალიდან:

slope(უჯრა) = (უჯრის სიმაღლე − თვალის სიმაღლე) / მანძილი

სამიზნე ჩანს მაშინ და მხოლოდ მაშინ, თუ მისი დახრილობა ყველა შუალედური უჯრის დახრილობაზე არანაკლებია. ასე ვინახავთ მხოლოდ ერთ რიცხვს, ზრდადი მაქსიმუმის სახით (max_slope), და ყოველ უჯრაზე ერთი შედარება გვაქვს.

უჯრებს ხაზზე ბრეზენჰემის ალგორითმი (Bresenham) ირჩევს: ის რასტრზე სწორი ხაზის უახლოეს უჯრებს აბრუნებს მხოლოდ მთელი რიცხვებით.

🎮 სცადეთ თქვენც

დააჭირეთ რუკას დამკვირვებლის გადასაადგილებლად. კურსორის ქვეშ ხედვის ხაზი გაივლება და ქვემოთ გამოჩნდება რელიეფის პროფილი.

დაკვირვება

  • დამკვირვებლის სიმაღლის ზრდა (კოშკური, ანძა) ხედვის არეს იმატებს, თუმცა ერთნაირი სისწრაფით არა: პირველი ათი მეტრი უფრო ეფექტურია, ვიდრე მომდევნო ორმოცი (იხ. დავალება 1).
  • სამიზნის სიმაღლე ეხმარება: შენობა ან ანძა მაშინაც ჩანს, თუ მისი ძირი ჩრდილშია.

🐍 რეალიზაცია სუფთა Python-ით

view.py
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
import math


def bresenham(r0, c0, r1, c1):
    """უჯრები (r0,c0)-დან (r1,c1)-მდე სწორ ხაზზე, ბოლოების ჩათვლით."""
    cells = []
    dr, dc = abs(r1 - r0), abs(c1 - c0)
    sr = 1 if r1 > r0 else -1
    sc = 1 if c1 > c0 else -1
    err = dc - dr
    r, c = r0, c0
    while True:
        cells.append((r, c))
        if (r, c) == (r1, c1):
            return cells
        e2 = 2 * err
        if e2 > -dr:
            err -= dr
            c += sc
        if e2 < dc:
            err += dc
            r += sr


def line_of_sight(dem, obs, target, h_obs=1.7, h_target=0.0, cell=30.0):
    """ხედავს თუ არა დამკვირვებელი obs=(r,c) სამიზნეს target=(r,c).
    cell = უჯრის ზომა მეტრებში, სიმაღლეები მეტრებში."""
    (r0, c0), (r1, c1) = obs, target
    z0 = dem[r0][c0] + h_obs
    max_slope = -math.inf
    path = bresenham(r0, c0, r1, c1)[1:]
    for r, c in path:
        d = math.hypot(r - r0, c - c0) * cell
        if (r, c) == target:
            return (dem[r][c] + h_target - z0) / d >= max_slope
        max_slope = max(max_slope, (dem[r][c] - z0) / d)   # ხაზი ამ უჯრის ზედაპირს ეხება
    return True                                            # obs == target


def viewshed(dem, obs, h_obs=1.7, h_target=0.0, cell=30.0):
    """ლოგიკური მატრიცა: True = უჯრა ჩანს."""
    rows, cols = len(dem), len(dem[0])
    return [[line_of_sight(dem, obs, (r, c), h_obs, h_target, cell) for c in range(cols)] for r in range(rows)]

გამოყენება

import math
from view import *


def terr(r, c):
    return 120 + 60 * math.sin(c / 9) * math.cos(r / 7) + 35 * math.sin((r + c) / 13) + 25 * math.cos(c / 4.5) * math.sin(r / 5)


dem = [[terr(r, c) for c in range(60)] for r in range(40)]     # 40 x 60, უჯრა 30 მ
vs = viewshed(dem, (20, 30), h_obs=10)
print('ხილულია:', sum(map(sum, vs)), 'უჯრა 2400-დან')
print(bresenham(0, 0, 3, 7))

ხილულია: 1091 უჯრა 2400-დან
[(0, 0), (0, 1), (1, 2), (1, 3), (2, 4), (2, 5), (3, 6), (3, 7)]

შემოწმება

სამი მარტივი შემთხვევა ისე იქცევა, როგორც უნდა: სრულიად ბრტყელ რელიეფზე ყველა უჯრა ჩანს; ერთ კედელზე (სიმაღლე 100) თავად კედელი ჩანს, მის უკან მყოფი უჯრები კი არა, ხოლო კედლამდე ყველაფერი ჩანს. JavaScript-ის ვერსია, რომელიც დემოშია, იგივე რიცხვებს იძლევა (1000, 1091, 1385 უჯრა 1.7, 10 და 50 მ-ზე).

სირთულე

ყოველი უჯრისთვის ხაზი O(L) ბიჯია (L რასტრის სიგრძეა), ამიტომ მთელი ხედვის არე O(N·L). ეს მარტივი „ყველა ხაზი ცალ-ცალკე" ვერსიაა. R2/R3 და sweep-line ალგორითმები (მაგ. GRASS-ში) ხედვის არეს გაცილებით სწრაფად იღებენ, ზოგს ინფორმაციას ხაზებს შორის ინახავენ.

📚 იგივე ბიბლიოთეკით

gdal_viewshed -ox 44.80 -oy 41.72 -oz 10 -tz 0 -md 20000 dem.tif viewshed.tif

-oz დამკვირვებლის სიმაღლეა, -tz სამიზნის, -md მაქსიმალური მანძილი. ნაგულისხმევად GDAL ითვალისწინებს დედამიწის სიმრუდესაც (-cc კოეფიციენტი).

import processing

processing.run('gdal:viewshed', {
    'INPUT': 'dem.tif', 'BAND': 1,
    'OBSERVER': '44.80,41.72 [EPSG:4326]',
    'OBSERVER_HEIGHT': 10, 'TARGET_HEIGHT': 0,
    'MAX_DISTANCE': 20000, 'OUTPUT': 'viewshed.tif'})

რა არის შემოწმებული

ბიბლიოთეკების ფრაგმენტები არ გამიშვია ამ გარემოში (GDAL და QGIS არ არის დაყენებული). პარამეტრების სახელები დოკუმენტაციასთან შეამოწმეთ, განსაკუთრებით GDAL-ის ვერსიაზე, რადგან gdal_viewshed ახალი ბრძანებაა.

✍️ დავალებები

დავალება 1: სიმაღლის გავლენა

გამოთვალეთ ხილული უჯრების წილი (20, 30)-დან დამკვირვებლის სიმაღლეებზე 0, 2, 10, 50 მეტრი.

ამოხსნა 

total = 60 * 40
for h in (0, 2, 10, 50):
    n = sum(map(sum, viewshed(dem, (20, 30), h)))
    print(f'სიმაღლე {h:2} მ: {n:4} უჯრა ({100 * n / total:.1f}%)')

სიმაღლე  0 მ:  958 უჯრა (39.9%)
სიმაღლე  2 მ: 1001 უჯრა (41.7%)
სიმაღლე 10 მ: 1091 უჯრა (45.5%)
სიმაღლე 50 მ: 1385 უჯრა (57.7%)
ზრდა არაწრფივია: 0-დან 10 მეტრამდე +5.6 პროცენტული პუნქტი, 10-დან 50-მდე +12.2, თუმცა ყოველი მეტრი ნაკლებს იძლევა (ეფექტი სიმაღლის მეტრზე: ~0.56 პ.პ. პირველ 10 მ-ზე, ~0.31 მომდევნო 40 მ-ზე). ამიტომ „კიდევ ერთი მეტრი" გვიან ნაკლებად ღირს.

დავალება 2: არის თუ არა ხედვა სიმეტრიული?

თუ A ხედავს B-ს, ხედავს თუ არა B A-ს? შეამოწმეთ 2 000 შემთხვევით წყვილზე (ორივე სიმაღლე 0).

ამოხსნა 

import random
random.seed(1)
asym = 0
N = 2000
for _ in range(N):
    a = (random.randrange(40), random.randrange(60))
    b = (random.randrange(40), random.randrange(60))
    if line_of_sight(dem, a, b, 0, 0) != line_of_sight(dem, b, a, 0, 0):
        asym += 1
print('ასიმეტრიული წყვილი:', asym, 'of', N)

ასიმეტრიული წყვილი: 28 of 2000
თეორიულად ხედვა სიმეტრიულია (ერთსა და იმავე ხაზზე ორივე მიმართულებით), მაგრამ რასტრზე არა: A-დან B-მდე და B-დან A-მდე ბრეზენჰემი სხვადასხვა უჯრებს ირჩევს, ამიტომ 1.4% წყვილში პასუხი განსხვავდება. თუ სიმეტრია გჭირდებათ (მაგ. ორმხრივი რადიოკავშირი), ორივე მიმართულება შეამოწმეთ და „და" ოპერაცია გამოიყენეთ. დამკვირვებლისა და სამიზნის სხვადასხვა სიმაღლეც ასიმეტრიას იძლევა.

დავალება 3: დედამიწის სიმრუდე

დაამატეთ line_of_sight-ს ვერსია, რომელიც მანძილზე d სიმაღლეს ამცირებს (1 − k)·d² / (2R)-ით (R = 6 371 km, რეფრაქციის k = 0.13). შეადარეთ ხილული უჯრები უჯრის ზომებზე 30, 300 და 1000 მ (იგივე რელიეფის სიმაღლეებით).

ამოხსნა 

R_EARTH = 6371000.0
K_REFR = 0.13


def los_curved(dem, obs, target, h_obs, h_target=0.0, cell=30.0):
    (r0, c0), (r1, c1) = obs, target
    z0 = dem[r0][c0] + h_obs
    drop = lambda d: (1 - K_REFR) * d * d / (2 * R_EARTH)
    max_slope = -math.inf
    for r, c in bresenham(r0, c0, r1, c1)[1:]:
        d = math.hypot(r - r0, c - c0) * cell
        if (r, c) == target:
            return (dem[r][c] + h_target - drop(d) - z0) / d >= max_slope
        max_slope = max(max_slope, (dem[r][c] - drop(d) - z0) / d)
    return True


obs = (20, 30)
for cell in (30, 300, 1000):
    flat = sum(sum(line_of_sight(dem, obs, (r, c), 10, 0, cell) for c in range(60)) for r in range(40))
    curv = sum(sum(los_curved(dem, obs, (r, c), 10, 0, cell) for c in range(60)) for r in range(40))
    print(f'უჯრა {cell:4} მ: სიმრუდის გარეშე {flat}, სიმრუდით {curv}')

უჯრა   30 მ: სიმრუდის გარეშე 1091, სიმრუდით 1090
უჯრა  300 მ: სიმრუდის გარეშე 1091, სიმრუდით 1076
უჯრა 1000 მ: სიმრუდის გარეშე 1091, სიმრუდით 870
ორი საინტერესო დასკვნა. ბრტყელ დედამიწაზე ხილვადობა მასშტაბზე არ არის დამოკიდებული (ყველა მანძილი ერთნაირად იზრდება, ამიტომ დახრილობათა შედარება არ იცვლება). სიმრუდე კი მასშტაბს აფასებს: 30 მეტრიან უჯრაზე (1.8 კმ-იანი რუკა) მისი გავლენა უმნიშვნელოა, 1 კმ-იანზე (60 კმ-იანი რუკა) კი ხედვის არე ~20%-ით მცირდება. ამიტომ დიდ მანძილებზე სიმრუდე აუცილებლად უნდა ჩაითვალოს.

დავალება 4: სად დავდგათ კოშკურები?

აირჩიეთ 4 წერტილი (ბადე ყოველ 5 უჯრაზე, სიმაღლე 10 მ) ისე, რომ ერთად რაც შეიძლება მეტი ტერიტორია დაფარონ. გამოიყენეთ ხარბი ალგორითმი: ყოველ ბიჯზე ის წერტილი, რომელიც ჯერ დაუფარავს ყველაზე მეტს მატებს.

ამოხსნა 

def greedy(dem, k, h=10, step=5):
    cands = [(r, c) for r in range(0, 40, step) for c in range(0, 60, step)]
    vs = {p: {(r, c) for r in range(40) for c in range(60) if line_of_sight(dem, p, (r, c), h)} for p in cands}
    chosen, covered = [], set()
    for _ in range(k):
        best = max(cands, key=lambda p: len(vs[p] - covered))
        chosen.append(best)
        covered |= vs[best]
        print(f'  + {best}: ჯამური დაფარვა {100 * len(covered) / total:.1f}%')
    return chosen


greedy(dem, 4)

  + (20, 40): ჯამური დაფარვა 64.9%
  + (35, 55): ჯამური დაფარვა 84.1%
  + (0, 15): ჯამური დაფარვა 93.7%
  + (10, 55): ჯამური დაფარვა 96.5%
ეს maximum coverage ამოცანაა: ზუსტი ამოხსნა NP-რთულია, ხარბი კი გარანტირებულად ოპტიმუმის ≥ 63%-ს იძლევა და პრაქტიკაში ხშირად ბევრად უკეთესია. პირველი კოშკურა საუკეთესო სივრცეს ფარავს, შემდეგი თითოეული ნაკლებს: ეს კლებადი სარგებლის კლასიკური სურათია.