Visualize Viewsheds and Coverage Maps Using Terrain

lat1 = 34.076389;
lon1 = -118.468944;
towerHeight1 = 17.1;
text1 = "Tower 1";

lat2 = 34.386389;
lon2 = -118.329778;
towerHeight2 = 30.5;
text2 = "Tower 2";

figure
basemap = "satellite";
geobasemap(basemap)

geoplot([lat1 lat2],[lon1 lon2], ...
    LineStyle="none",Marker="o",MarkerSize=10, ...
    MarkerEdgeColor="k",MarkerFaceColor="c")

gx = gca;
zoomLevel = floor(gx.ZoomLevel);
gx.ZoomLevel = zoomLevel;        

[latlim,lonlim] = geolimits(gx);

txtDelta = 0.005;
text(lat1 + txtDelta,lon1 + txtDelta,text1,Color="w")
text(lat2 + txtDelta,lon2 + txtDelta,text2,Color="w")
title("Satellite Imagery of San Fernando Valley Region")
subtitle("With Placement of Cell Towers")

layer = wmsfind("mathworks",SearchField="serverurl");
layer = refine(layer,"elevation");

[Z,RZ] = wmsread(layer,Latlim=latlim,Lonlim=lonlim,ImageFormat="image/bil");

Z = double(Z);

height1 = geointerp(Z,RZ,lat1,lon1,"nearest")

height2 = geointerp(Z,RZ,lat2,lon2,"nearest")

figure
usamap(Z,RZ)
geoshow(Z,RZ,DisplayType="texturemap")
demcmap(Z,256)

geoshow(lat1,lon1,DisplayType="point",ZData=height1, ...
    MarkerEdgeColor="k",MarkerFaceColor="c",MarkerSize=10,Marker="o")
geoshow(lat2,lon2,DisplayType="point",ZData=height2, ...
     MarkerEdgeColor="k",MarkerFaceColor="c",MarkerSize=10,Marker="o")

textm(lat1 + txtDelta,lon1 + txtDelta,text1)
textm(lat2 + txtDelta,lon2 + txtDelta,text2)
title("Elevation of San Fernando Valley Region")

cb = colorbar;
cb.Label.String = "Elevation in meters";

[A,RA,attrib] = readBasemapImage(basemap,latlim,lonlim,zoomLevel);

proj = RA.ProjectedCRS;
[x1,y1] = projfwd(proj,lat1,lon1);
[x2,y2] = projfwd(proj,lat2,lon2);

figure
axis off
hold on

mapshow(A,RA)
plot([x1 x2],[y1 y2], ...
    LineStyle="none",Marker="o",MarkerSize=10, ...
    MarkerEdgeColor="k",MarkerFaceColor="c")

txtDeltaXY = 500;
text(x1 + txtDeltaXY,y1 + txtDeltaXY,text1,Color="w")
text(x2 + txtDeltaXY,y2 + txtDeltaXY,text2,Color="w")
title("Satellite Basemap Imagery of San Fernando Valley Region")
subtitle("Attribution: " + attrib)

[latlim,lonlim] = projinv(proj,RA.XWorldLimits,RA.YWorldLimits);
szA = RA.RasterSize;

[Z,RZ] = wmsread(layer,Latlim=latlim,Lonlim=lonlim, ...
    ImageHeight=szA(1),ImageWidth=szA(2),ImageFormat="image/bil");
Z = double(Z);

[vis1,visR1] = viewshed(Z,RZ,lat1,lon1,towerHeight1);
[vis2,visR2] = viewshed(Z,RZ,lat2,lon2,towerHeight2);

visCombined = vis1 | vis2;
vis = visCombined;
vis = double(vis);
vis(vis == 0) = NaN;
visR = visR2;

[latz,lonz] = geographicGrid(RZ);
[xz,yz] = projfwd(proj,latz,lonz);

cmapt = [0 0 0];
covcolor = [0.1034  0.8960  0.7150];
cmapt(2,:,:) = covcolor;

figure
clim([0 1])
colormap(cmapt)
tiledlayout(2,3,Padding="tight",TileSpacing="tight")

nexttile
geoshow(vis1,visR1,EdgeColor="none",DisplayType="surface")
axis off
title("Viewshed")
subtitle("Tower 1")

nexttile(4)
geoshow(vis2,visR2,EdgeColor="none",DisplayType="surface")
axis off
title("Viewshed")
subtitle("Tower 2")

nexttile([2,2])
mapshow(A,RA)
hold on
mapshow(xz,yz,vis,EdgeColor="none",DisplayType="surface")
axis off

zdelta = 150;
z1 = height1 + towerHeight1 + zdelta;
z2 = height2 + towerHeight2 + zdelta;

text(x1,y1,z1,text1,Color="k",EdgeColor="k",BackgroundColor="w")
text(x2,y2,z2,text2,Color="k",EdgeColor="k",BackgroundColor="w")
title("Combined Viewshed Over Basemap Image")
subtitle("Attribution: " + attrib)

[row,col] = find(visCombined);
sz = size(visCombined);
ind = sub2ind(sz,row,col);
vis(ind) = Z(ind);

figure
colormap(covcolor)

hold on
axis off 
axis tight

mapshow(xz,yz,Z,DisplayType="surface",CData=A)

mapshow(xz,yz,vis,EdgeColor="none",DisplayType="surface")

text(x1,y1,z1,text1,Color="k",EdgeColor="k",BackgroundColor="w")
text(x2,y2,z2,text2,Color="k",EdgeColor="k",BackgroundColor="w")
title("Satellite Basemap Imagery with Terrain and Viewshed")
subtitle("Attribution: " + attrib)

daspect([1 1 0.25])
view(3)

txTower1 = txsite(Name=text1,Latitude=lat1,Longitude=lon1);
txTower2 = txsite(Name=text2,Latitude=lat2,Longitude=lon2);
txsites = [txTower1; txTower2];
pd = coverage(txsites);

dataVariableName = "Power";
maxrange = [30000 30000];
[lats,lons,data,levels] = propagationDataGrid(pd,dataVariableName,maxrange,txsites);
[covlatlim,covlonlim] = geoquadline(lats,lons);

cmap = turbo(length(levels));
climLevels = [min(levels) max(levels)+10];

figure
colormap(cmap)
clim(climLevels)
hold on

usamap(covlatlim,covlonlim)
geoshow(lats,lons,data,DisplayType="surface")
geoshow(latlim([1 2 2 1 1]),lonlim([1,1,2,2,1]),Color="k")

textm(lat1,lon1,text1,Color="k",FontSize=12)
textm(lat2,lon2,text2,Color="k",FontSize=12)
title("Coverage Map with Boundary of Basemap Image")

cb = colorbar;
cb.Label.String = "Power (dBm)";
cb.Ticks = levels;

[cx,cy] = projfwd(proj,lats,lons);

figure
colormap(cmap)
clim(climLevels)
axis off
hold on

mapshow(A,RA)
contourf(cx,cy,data,LevelList=levels)

xlim(RA.XWorldLimits)
ylim(RA.YWorldLimits)

text(x1,y1,text1,Color="w") 
text(x2,y2,text2,Color="w")
title("Contoured Coverage Map Over Basemap Image")
subtitle("Attribution: " + attrib)

cb = colorbar;
cb.Label.String = "Power (dBm)";
cb.Ticks = levels;

[vlat,vlon] = geographicGrid(visR);
[vx,vy] = projfwd(proj,vlat,vlon);

figure
ax1 = axes;
hold(ax1,"on")
axis(ax1,"off")
clim(ax1,climLevels)
colormap(ax1,cmap)

mapshow(ax1,A,RA)
contourf(ax1,cx,cy,data,LevelList=levels);

xlim(ax1,RA.XWorldLimits)
ylim(ax1,RA.YWorldLimits)
title(ax1,"Contoured Coverage Map and Viewshed on Basemap Image")
subtitle("Attribution: " + attrib)

cb1 = colorbar(ax1);
cb1.Label.String = "Power (dBm)";
cb1.Ticks = levels;

ax2 = axes(Visible="off");
hold(ax2,"on")
colormap(ax2,covcolor)

[~,visContourMap] = contourf(ax2,vx,vy,vis,LevelList=0.5:1:1.5);

cb2 = colorbar(ax2);
cb2.Label.String = "Viewshed";
cb2.Ticks = [];

ax1.Position = [ax1.Position(1) ax1.Position(2)-0.04 ax1.Position(3) ax1.Position(4)];
ax2.Position = ax1.Position;
cb1.Position = [cb1.Position(1) cb1.Position(2) cb1.Position(3) cb1.Position(4)*0.8];
cb2.Position = [cb2.Position(1) 0.76 cb2.Position(3) 0.1];

faceAlpha = 0.3;
visContourMap.FaceAlpha = faceAlpha;

GT = contourDataGrid(lats,lons,data,levels);
disp(GT)

figure
geobasemap satellite
colormap(cmap)
clim(climLevels)
hold on

hp = geoplot(GT,ColorVariable="Power (dBm)",EdgeColor="k",FaceAlpha=1);
boundary = geoplot(latlim([1 2 2 1 1]),lonlim([1,1,2,2,1]),Color="k",LineWidth=2);

text(lat1,lon1,text1,Color="w") 
text(lat2,lon2,text2,Color="w")
title("Contoured Coverage Map with Boundary of Basemap Image")

cb = colorbar;
cb.Label.String = "Power (dBm)";
cb.Ticks = levels;

geolimits(latlim,lonlim)
delete(boundary)
title("Contoured Coverage Map at Boundary of Basemap Image")

faceAlpha =0.3;
hp.FaceAlpha = faceAlpha;

figure
colormap(cmap)
clim(climLevels)
hold on

bar(GT.("Power (dBm)"),GT.("Area (sq km)"),FaceColor="flat",CData=cmap);
ylabel("Area in square kilometers")
xlabel("Power (dBm)")

cb = colorbar;
cb.Label.String = "Power (dBm)";
cb.Ticks = levels;

locationLats = [34.2107 34.1493];
locationLons = [-118.4545 -118.1744];

[tf,powerLevels] = incoverage(locationLats,locationLons,GT);
disp(tf)

disp(powerLevels)

figure
geobasemap topographic
hold on
geotickformat dd
geolimits(latlim,lonlim)

geoplot(locationLats(tf),locationLons(tf),"o", ...
    MarkerEdgeColor="c",MarkerFaceColor="c",MarkerSize=10)
geoplot(locationLats(~tf),locationLons(~tf),"o", ...
    MarkerEdgeColor="k",MarkerFaceColor="k",MarkerSize=10)

figure
geobasemap topographic
geolimits(latlim,lonlim)
geotickformat dd
hold on
colormap(cmap)
clim(climLevels)

contourPlot = geoplot(GT,ColorVariable="Power (dBm)",EdgeColor="k",FaceAlpha=1);

cb = colorbar;
cb.Label.String = "Power (dBm)";
cb.Ticks = levels(1:end);

interactivelySelectPoints =false;
if interactivelySelectPoints
    numberOfPoints = 4; %#ok<*UNRCH> 
    title("Select " + string(numberOfPoints) + " Points")
    pts = ginput(numberOfPoints);
else
    pts = [ ...
        34.0732 -118.4652
        34.1880 -118.2000
        34.2200 -118.6172
        34.3849 -118.5472];
end

for k = 1:size(pts,1)
    latpoint = pts(k,1);
    lonpoint = pts(k,2);
    [tf,powerLevel] = incoverage(latpoint,lonpoint,GT);

    if ~tf
        color = [0 0 0];
    else
        index = GT.("Power (dBm)") == powerLevel;
        color = cmap(index,:,:);
    end

    h = geoplot(latpoint,lonpoint,"o",MarkerEdgeColor=color,MarkerFaceColor=color,MarkerSize=10);
    dt = datatip(h);
    dtrow = dataTipTextRow("Power",powerLevel);
    h.DataTipTemplate.DataTipRows(end+1) = dtrow;
end
title("Selected Points")

showContourPlot = false;
contourPlot.Visible = showContourPlot;

T = geocoverage(txsites);
disp(T)

basemap = "topographic";
geocoverage(txsites(1),basemap)

names = ["Clarendon St.","Cambridge St.","Albany St."];
bostonLat = [42.348722 42.361222 42.338444];
bostonLon = [-71.075889 -71.069778 -71.065611];
bostonH = [260 30 23];
freq = [852.637e6 862.012e6 862.012e6];
txs = txsite(Name=names,Latitude=bostonLat,Longitude=bostonLon, ...
             AntennaHeight=bostonH,TransmitterFrequency=freq);
geocoverage(txs,basemap)

gx = gca;
gx.ZoomLevel = 12;

function [clats,clons,cdata,levels] = propagationDataGrid(pd,dataVariableName,maxrange,txsites)
% Extract propagation data grid from propagationData object

    % Create grids of latitude and longitude values
    imageSize = [512 512];
    lats = pd.Data.Latitude;
    lons = pd.Data.Longitude;
    data = pd.Data.(dataVariableName);

    [lonmin,lonmax] = bounds(lons);
    [latmin,latmax] = bounds(lats);
    lonsv = linspace(lonmin,lonmax,imageSize(2));
    latsv = linspace(latmin,latmax,imageSize(1));
    [clons,clats] = meshgrid(lonsv,latsv);

    % Determine which latitude and longitude values are within range of the
    % transmitters
    txlat = [txsites.Latitude]';
    txlon = [txsites.Longitude]';

    numTx = size(txlat,1);
    gc = nan(numel(clons),numTx);
    for txInd = 1:numTx
        gc(:,txInd) = distance(txlat(txInd),txlon(txInd),clats(:),clons(:),wgs84Ellipsoid);
    end
    isInRange = any(gc <= maxrange,2);

    % Create grid of power values
    cdata = nan(imageSize);
    cdata(isInRange) = interp(pd,clats(isInRange),clons(isInRange), ...
        DataVariableName=dataVariableName);

    % Determine power levels
    [bmin,bmax] = bounds(data);
    levels = max(bmin,-120):10:bmax;  
    datalevels = sort(levels);
    maxBin = max(max(datalevels(:)),max(cdata(:))) + 1; % Need max bin edge to include all values
    bins = [datalevels(:); maxBin];
    cdata = discretize(cdata,bins,datalevels);
end

function GT = contourDataGrid(latd,lond,data,levels)
% Contour data grid created from propagationDataGrid helper function

    % Replace NaN values in power grid with a large negative number
    data(isnan(data)) = min(levels) - 1000;

    % Project the coordinates using an equal-area projection
    proj = projcrs(102008,"Authority","ESRI");
    [xd,yd] = projfwd(proj,latd,lond);

    % Contour the projected data
    fig = figure("Visible","off");
    obj = onCleanup(@()close(fig));
    c = contourf(xd,yd,data,LevelList=levels);

    % Initialize variables
    x = c(1,:);
    y = c(2,:);
    n = length(y);
    k = 1;
    index = 1;
    levels = zeros(n,1);
    latc = cell(n,1);
    lonc = cell(n,1);

    % Process the coordinates of the projected contours
    while k < n
        level = x(k);
        numVertices = y(k);
        yk = y(k+1:k+numVertices);
        xk = x(k+1:k+numVertices);
        k = k + numVertices + 1;

        [first,last] = findNanDelimitedParts(xk);
        jindex = 0;
        jy = {};
        jx = {};

        for j = 1:numel(first)
            % Process the j-th part of the k-th level
            s = first(j);
            e = last(j);
            cx = xk(s:e);
            cy = yk(s:e);
            if cx(end) ~= cx(1) && cy(end) ~= cy(1)
                cy(end+1) = cy(1); %#ok<*AGROW>
                cx(end+1) = cx(1);
            end

            if j == 1 && ~ispolycw(cx,cy)
                % First region must always be clockwise
                [cx,cy] = poly2cw(cx,cy);
            end

            % Filter regions made of collinear points or fewer than 3 points
            if ~(length(cx) < 4 || all(diff(cx) == 0) || all(diff(cy) == 0))
                jindex = jindex + 1;
                jy{jindex,1} = cy(:)';
                jx{jindex,1} = cx(:)';
            end
            
        end

        % Unproject the coordinates of the processed contours
        [jx,jy] = polyjoin(jx,jy);
        if length(jy) > 2 && length(jx) > 2
            [jlat,jlon] = projinv(proj,jx,jy);
            latc{index,1} = jlat(:)';
            lonc{index,1} = jlon(:)';
            levels(index,1) = level;
            index = index + 1;
        end
    end

    % Create contour shapes from the unprojected coordinates. Include the
    % contour shapes, the areas of the shapes, the perimeters of the
    % shapes, and the power levels in a geospatial table.
    latc = latc(1:index-1);
    lonc = lonc(1:index-1);
    Shape = geopolyshape(latc,lonc);

    Area = area(Shape);

    Perimeter = perimeter(Shape);

    levels = levels(1:index-1);
    Levels = levels;

    allPartsGT = table(Shape,Area,Perimeter,Levels);  

    % Condense the geospatial table into a new geospatial table with one
    % row per power level.
    GT = table.empty;
    levels = unique(allPartsGT.Levels);
    for k = 1:length(levels)
        gtLevel = allPartsGT(allPartsGT.Levels == levels(k),:);
        tLevel = geotable2table(gtLevel,["Latitude","Longitude"]);
        [lon,lat] = polyjoin(tLevel.Longitude,tLevel.Latitude);
        Shape = geopolyshape(lat,lon);
        Levels = levels(k);
        Area = area(Shape);
        Perimeter = perimeter(Shape);
        GT(end+1,:) = table(Shape,Area,Perimeter,Levels);
    end

    maxLevelDiff = max(abs(diff(GT.Levels)));
    LevelRange = [GT.Levels GT.Levels + maxLevelDiff];
    GT.LevelRange = LevelRange;

    % Clean up the geospatial table
    GT.Area = GT.Area/1e6;
    GT = sortrows(GT,"Area","descend");
    GT.Properties.VariableNames = ...
        ["Shape","Area (sq km)","Perimeter (km)","Power (dBm)","Power Range (dBm)"]; 
end

function [in,powerLevels] = incoverage(lat,lon,T)
% Query points in coverage

    % Determine whether points are within coverage
    tf = false(length(T.Shape),length(lat));
    for k = 1:length(T.Shape)
        tf(k,:) = isinterior(T.Shape(k),geopointshape(lat,lon));
    end

    % Determine the power levels for the points
    in = false(length(lat),1);
    powerLevels = -inf(length(lat),1);
    for k = 1:length(in)
        v = tf(:,k);
        in(k) = any(v);
        if in(k)
            powerLevels(k) = max(T{v,"Power (dBm)"});
        end
    end
end

function varargout = geocoverage(txsites,basemap)
% Display or calculate contoured coverage map

    % Calculate and contour the coverage map
    pd = coverage(txsites);
    dataVariableName = "Power";
    maxrange = 30000*ones(1,length(txsites));
    [lats,lons,data,levels] = propagationDataGrid(pd,dataVariableName,maxrange,txsites);
    GT = contourDataGrid(lats,lons,data,levels);

    % If you do not specify an output argument, display the coverage map
    if nargout == 0
        if nargin < 2
            basemap = "satellite";
        end

        cmap = turbo(length(levels));
        figure
        geobasemap(basemap)

        geoplot(GT,ColorVariable="Power (dBm)",EdgeColor="k",FaceAlpha=0.5)
        colormap(cmap)
        clim([min(levels) max(levels)+10])

        h = colorbar;
        h.Label.String = "Power (dBm)";
        h.Ticks = levels;

        for k = 1:length(txsites)
            latt = txsites(k).Latitude;
            lont = txsites(k).Longitude;
            text(latt,lont,txsites(k).Name,Color="w")
        end

        title("Coverage Contour Map")
    % If you specify an output argument, return a geospatial table
    else
        varargout{1} = GT;
    end
end

function [first,last] = findNanDelimitedParts(x)
% Find indices of the first and last elements of each part in x. 
% x can contain runs of multiple NaNs, and x can start or end with 
% one or more NaNs.

    n = isnan(x(:));

    firstOrPrecededByNaN = [true; n(1:end-1)];
    first = find(~n & firstOrPrecededByNaN);

    lastOrFollowedByNaN = [n(2:end); true];
    last = find(~n & lastOrFollowedByNaN);
end