Coverage Maps for Satellite Constellation

This example uses:

This example shows how to display coverage maps for a region of interest using a satellite scenario and 2-D maps. A coverage map for satellites includes received power over an area for a downlink to ground-based receivers within the area.

By combining the mapping capabilities in Mapping Toolbox™ with the satellite scenario capabilities in Satellite Communications Toolbox, you can create contoured coverage maps and use the contour data for visualization and analysis. This example shows you how to:

Compute coverage map data using a satellite scenario containing a grid of ground station receivers.
View coverage maps using different map displays, including axesm-based maps and map axes.
Compute and view coverage maps for other times of interest.

Create and Visualize Scenario

Create a satellite scenario and a viewer. Specify the start date and duration of 1 hour for the scenario.

startTime = datetime(2023,2,21,18,0,0);
stopTime = startTime + hours(1);
sampleTime = 60; % seconds
sc = satelliteScenario(startTime,stopTime,sampleTime);
viewer = satelliteScenarioViewer(sc,ShowDetails=false);

Add Satellite Constellation

The Iridium NEXT satellite network, launched between 2018 and 2019 [1], contains 66 active LEO satellites with:

Six orbital planes with an approximate inclination of 86.6 degrees and difference of RAAN of approximately 30 degrees between planes [1].
11 satellites per orbital plane with an approximate difference in true anomaly of 32.7 degrees between satellites [1].

Add the active Iridium NEXT satellites to the scenario. Create the orbital elements for the satellites in the Iridium network and create the satellites.

numSatellitesPerOrbitalPlane = 11;
numOrbits = 6;

orbitIdx = repelem(1:numOrbits,1,numSatellitesPerOrbitalPlane);
planeIdx = repmat(1:numSatellitesPerOrbitalPlane,1,numOrbits);

RAAN = 180*(orbitIdx-1)/numOrbits;
trueanomaly = 360*(planeIdx-1 + 0.5*(mod(orbitIdx,2)-1))/numSatellitesPerOrbitalPlane;
semimajoraxis = repmat((6371 + 780)*1e3,size(RAAN)); % meters
inclination = repmat(86.4,size(RAAN)); % degrees
eccentricity = zeros(size(RAAN)); % degrees
argofperiapsis = zeros(size(RAAN)); % degrees

iridiumSatellites = satellite(sc,...
    semimajoraxis,eccentricity,inclination,RAAN,argofperiapsis,trueanomaly,...
    Name="Iridium " + string(1:66)');

If you have a license for Aerospace Toolbox™, you can also create the Iridium constellation using the walkerStar (Aerospace Toolbox) function.

Add Grid of Ground Stations Covering Australia

Create the coverage map data by computing signal strength for a grid of ground station receivers covering Australia. Create ground station locations by creating a grid of latitude-longitude coordinates and selecting the coordinates within a buffered region corresponding to mainland Australia.

Specify approximate latitude and longitude limits for Australia and a spacing in degrees. The spacing determines the distance between the ground stations. A lower spacing value results in improved coverage map resolution but also increases computation time.

latlim = [-40 -9];
lonlim = [110 160];
spacingInLatLon = 1; % degrees

Create a projected coordinate reference system (CRS) that is appropriate for Australia. A projected CRS provides information that assigns Cartesian x and y map coordinates to physical locations. Use the GDA2020 / GA LCC projected CRS appropriate for Australia that uses the Lambert Conformal Conic projection.

proj = projcrs(7845);

Calculate the projected X-Y map coordinates for the data grid. Use the projected X-Y map coordinates instead of latitude-longitude coordinates to reduce variation in the distance between grid locations.

spacingInXY = deg2km(spacingInLatLon)*1000; % meters
[xlimits,ylimits] = projfwd(proj,latlim,lonlim);
R = maprefpostings(xlimits,ylimits,spacingInXY,spacingInXY);
[X,Y] = worldGrid(R);

Transform the grid locations to latitude-longitude coordinates.

[gridlat,gridlon] = projinv(proj,X,Y);

Read a shapefile containing world land areas into the workspace as a geospatial table. The geospatial table represents the land areas using polygon shapes in geographic coordinates. Extract the polygon shape for Australia to define a buffered region.

landareas = readgeotable("landareas.shp");
australia = landareas(landareas.Name == "Australia",:);

Create a new polygon shape in geographic coordinates that includes a buffer region around the Australia land area. The buffer ensures the coverage map extents can cover the land area.

T = geotable2table(australia,["Latitude","Longitude"]);
[landlat,landlon] = polyjoin(T.Latitude,T.Longitude);
[landlat,landlon] = maptrimp(landlat,landlon,latlim,lonlim);

bufwidth = 1;
[landlatb,landlonb] = bufferm(landlat,landlon,bufwidth,"outPlusInterior");
australiab = geopolyshape(landlatb,landlonb);

Select the grid coordinates within the buffered land area region.

gridpts = geopointshape(gridlat,gridlon);
inregion = isinterior(australiab,gridpts);

gslat = gridlat(inregion);
gslon = gridlon(inregion);

Create ground stations for the coordinates.

gs = groundStation(sc,gslat,gslon);

Add Transmitters and Receivers

Enable the modeling of downlinks by adding transmitters to the satellites and adding receivers to the ground stations.

The Iridium constellation uses a 48-beam antenna. Select a Gaussian antenna or a custom 48-beam antenna (Phased Array System Toolbox™ required). The Gaussian antenna approximates the 48-beam antenna by using the Iridium half view angle of 62.7 degrees [1]. Orient the antenna to point towards the nadir direction.

fq = 1625e6; % Hz
txpower = 20; % dBW
antennaType = "Gaussian";
halfBeamWidth = 62.7; % degrees

Add satellite transmitters with Gaussian antennas.

if antennaType == "Gaussian"
    lambda = physconst('lightspeed')/fq; % meters
    dishD = (70*lambda)/(2*halfBeamWidth); % meters
    tx = transmitter(iridiumSatellites, ...
        Frequency=fq, ...
        Power=txpower); 
    gaussianAntenna(tx,DishDiameter=dishD);
end

Add satellite transmitters with custom 48-beam antenna.

if antennaType == "Custom 48-Beam"
    antenna = helperCustom48BeamAntenna(fq);
    tx = transmitter(iridiumSatellites, ...
        Frequency=fq, ...
        MountingAngles=[0,-90,0], ... % [yaw, pitch, roll] with -90 using Phased Array System Toolbox convention
        Power=txpower, ...
        Antenna=antenna);  
end

Add Ground Station Receivers

Add ground station receivers with isotropic antennas.

isotropic = arrayConfig(Size=[1 1]);
rx = receiver(gs,Antenna=isotropic);

Visualize the satellite transmitter antenna patterns. The visualization shows the nadir-pointing orientation of the antennas.

pattern(tx,Size=500000);

Compute Raster Coverage Map Data

Close satellite scenario viewer and compute coverage data corresponding to the scenario start time. The satcoverage helper function returns gridded coverage data. For each satellite, the function computes a field-of-view shape and calculates signal strength for each ground station receiver within the satellite field-of-view. The coverage data is the maximum signal strength received from any satellite.

delete(viewer)
maxsigstrength = satcoverage(gridpts,sc,startTime,inregion,halfBeamWidth);

Visualize Coverage on an `axesm`-Based Map

Plot the raster coverage map data on an axesm-based map corresponding to Australia. An axesm-based map is used because it supports plotting a mix of vector and raster data, whereas map axes supports plotting vector data only.

Define minimum and maximum power levels for the display.

minpowerlevel = -120; % dBm
maxpowerlevel = -106; % dBm

Create a world map for Australia and plot the raster coverage map data as a contour display.

figure
worldmap(latlim,lonlim)
mlabel south

colormap turbo
clim([minpowerlevel maxpowerlevel])
geoshow(gridlat,gridlon,maxsigstrength,DisplayType="contour",Fill="on")
geoshow(australia,FaceColor="none")

cBar = contourcbar;
title(cBar,"dBm");
title("Signal Strength at " + string(startTime) + " UTC")

sydlatlon = [-32.13 151.21]; % Sydney
mellatlot = [-36.19 144.96]; % Melbourne
brislatlon = [-26.53 153.03]; % Brisbane

textm(sydlatlon(1),sydlatlon(2),"Sydney")
textm(mellatlot(1),mellatlot(2),"Melbourne")
textm(brislatlon(1),brislatlon(2),"Brisbane")

Compute Coverage Map Contours

Contour the raster coverage map data to create a geospatial table of coverage map contours. Each coverage map contour is represented as a row in the geospatial table containing a polygon shape in geographic coordinates. You can use the coverage map contours for both visualization and analysis.

levels = linspace(minpowerlevel,maxpowerlevel,8);
GT = contourDataGrid(gridlat,gridlon,maxsigstrength,levels,proj);
GT = sortrows(GT,"Power (dBm)");
disp(GT)

       Shape        Area (sq km)    Power (dBm)    Power Range (dBm)
    ____________    ____________    ___________    _________________

    geopolyshape     8.5597e+06        -120          -120    -118   
    geopolyshape     8.2697e+06        -118          -118    -116   
    geopolyshape     5.6029e+06        -116          -116    -114   
    geopolyshape     3.7171e+06        -114          -114    -112   
    geopolyshape     2.2785e+06        -112          -112    -110   
    geopolyshape     1.1718e+06        -110          -110    -108   
    geopolyshape      3.706e+05        -108          -108    -106

Visualize Coverage on a Map Axes

Plot the coverage map contours on a map axes using the map projection corresponding to Australia. For more information, see Choose a 2-D Map Display in Mapping Toolbox.

figure
newmap(proj)
hold on

colormap turbo
clim([minpowerlevel maxpowerlevel])
geoplot(GT,ColorVariable="Power (dBm)",EdgeColor="none")
geoplot(australia,FaceColor="none")

cBar = colorbar;
title(cBar,"dBm");
title("Signal Strength at " + string(startTime) + " UTC")

text(sydlatlon(1),sydlatlon(2),"Sydney",HorizontalAlignment="center")
text(mellatlot(1),mellatlot(2),"Melbourne",HorizontalAlignment="center")
text(brislatlon(1),brislatlon(2),"Brisbane",HorizontalAlignment="center")

Figure contains an axes object with type mapaxes. The mapaxes object contains 5 objects of type polygon, text.

Compute and Visualize Coverage for a Different Time of Interest

Specify a different time of interest and compute and visualize coverage.

secondTOI = startTime + minutes(2); % 2 minutes after the start of the scenario
maxsigstrength = satcoverage(gridpts,sc,secondTOI,inregion,halfBeamWidth);

GT2 = contourDataGrid(gridlat,gridlon,maxsigstrength,levels,proj);
GT2 = sortrows(GT2,"Power (dBm)");

figure
newmap(proj)
hold on

colormap turbo
clim([minpowerlevel maxpowerlevel])
geoplot(GT2,ColorVariable="Power (dBm)",EdgeColor="none")
geoplot(australia,FaceColor="none")

cBar = colorbar;
title(cBar,"dBm");
title("Signal Strength at " + string(secondTOI) + " UTC")

text(sydlatlon(1),sydlatlon(2),"Sydney",HorizontalAlignment="center")
text(mellatlot(1),mellatlot(2),"Melbourne",HorizontalAlignment="center")
text(brislatlon(1),brislatlon(2),"Brisbane",HorizontalAlignment="center")

Figure contains an axes object with type mapaxes. The mapaxes object contains 5 objects of type polygon, text.

Compute Coverage Levels for Cities

Compute and display coverage levels for three cities in Australia.

covlevels1 = [coveragelevel(sydlatlon(1),sydlatlon(2),GT); ...
    coveragelevel(mellatlot(1),mellatlot(2),GT); ...
    coveragelevel(brislatlon(1),brislatlon(2),GT)];
covlevels2 = [coveragelevel(sydlatlon(1),sydlatlon(2),GT2); ...
    coveragelevel(mellatlot(1),mellatlot(2),GT2); ...
    coveragelevel(brislatlon(1),brislatlon(2),GT2)];

covlevels = table(covlevels1,covlevels2, ...
    RowNames=["Sydney" "Melbourne" "Brisbane"], ...
    VariableNames=["Signal Strength at T1 (dBm)" "Signal Strength T2 (dBm)"])

covlevels=3×2 table
                 Signal Strength at T1 (dBm)    Signal Strength T2 (dBm)
                 ___________________________    ________________________

    Sydney                  -108                          -114          
    Melbourne               -112                          -112          
    Brisbane                -112                          -118

Helper Functions

The satcoverage helper function returns coverage data. For each satellite, the function computes a field-of-view shape and calculates signal strength for each ground station receiver within the satellite field-of-view. The coverage data is the maximum signal strength received from any satellite.

function coveragedata = satcoverage(gridpts,sc,timeIn,inregion,halfBeamWidth)

    % Get satellites and ground station receivers
    sats = sc.Satellites;
    rxs = [sc.GroundStations.Receivers];

    % Compute the latitude, longitude, and altitude of all satellites at the input time
    lla = states(sats,timeIn,"CoordinateFrame","geographic");

    % Initialize coverage data
    coveragedata = NaN(size(gridpts));

    for satind = 1:numel(sats)
        % Create a geopolyshape for the satellite field-of-view
        fov = fieldOfViewShape(lla(:,1,satind),halfBeamWidth);

        % Find grid and rx locations which are within the field-of-view
        gridInFOV = isinterior(fov,gridpts);
        rxInFOV = gridInFOV(inregion);

        % Compute sigstrength at grid locations using temporary link objects
        gridsigstrength = NaN(size(gridpts));
        if any(rxInFOV)
            tx = sats(satind).Transmitters;
            lnks = link(tx,rxs(rxInFOV));
            rxsigstrength = sigstrength(lnks,timeIn)+30; % Convert from dBW to dBm
            gridsigstrength(inregion & gridInFOV) = rxsigstrength;
            delete(lnks)
        end

        % Update coverage data with maximum signal strength found
        coveragedata = max(gridsigstrength,coveragedata);
    end
end

The fieldOfViewShape helper function returns a geopolyshape object representing the field-of-view for a satellite position assuming nadir-pointing satellite and a spherical Earth model.

function satFOV = fieldOfViewShape(lla,halfViewAngle)

    % Find the Earth central angle using the beam view angle
    if isreal(acosd(sind(halfViewAngle)*(lla(3)+earthRadius)/earthRadius))
        % Case when Earth FOV is bigger than antenna FOV 
        earthCentralAngle = 90-acosd(sind(halfViewAngle)*(lla(3)+earthRadius)/earthRadius)-halfViewAngle;
    else
        % Case when antenna FOV is bigger than Earth FOV 
        earthCentralAngle = 90-halfViewAngle;
    end

    % Create a buffer zone centered at the position of the satellite with a buffer of width equaling the Earth central angle
    [latBuff,lonBuff] = bufferm(lla(1),lla(2),earthCentralAngle,"outPlusInterior",100);

    % Handle the buffer zone crossing over -180/180 degrees
    if sum(abs(lonBuff) == 180) > 2
        crossVal = find(abs(lonBuff)==180) + 1;
        lonBuff(crossVal(2):end) = lonBuff(crossVal(2):end) - 360 *sign(lonBuff(crossVal(2)));
    elseif sum(abs(lonBuff) == 180) == 2
        lonBuff = [lonBuff; lonBuff(end); lonBuff(1); lonBuff(1)];
        if latBuff(1) > 0
            latBuff = [latBuff; 90; 90; latBuff(1)];
        else
            latBuff = [latBuff; -90; -90; latBuff(1)];
        end
    end

    % Create geopolyshape from the resulting latitude and longitude buffer zone values
    satFOV = geopolyshape(latBuff,lonBuff);
end

The contourDataGrid helper function contours a data grid and returns the result as a geospatial table. Each row of the table corresponds to a power level and includes the contour shape, the area of the contour shape in square kilometers, the minimum power for the level, and the power range for the level.

function GT = contourDataGrid(latd,lond,data,levels,proj)

    % Pad each side of the grid to ensure no contours extend past the grid bounds
    lond = [2*lond(1,:)-lond(2,:); lond; 2*lond(end,:)-lond(end-1,:)];
    lond = [2*lond(:,1)-lond(:,2) lond 2*lond(:,end)-lond(:,end-1)];
    latd = [2*latd(1,:)-latd(2,:); latd; 2*latd(end,:)-latd(end-1,:)];
    latd = [2*latd(:,1)-latd(:,2) latd 2*latd(:,end)-latd(:,end-1)];
    data = [ones(size(data,1)+2,1)*NaN [ones(1,size(data,2))*NaN; data; ones(1,size(data,2))*NaN] ones(size(data,1)+2,1)*NaN];
    
    % 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
    [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, 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);
    
    levels = levels(1:index-1);
    Levels = levels;
    
    allPartsGT = table(Shape,Area,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);
        GT(end+1,:) = table(Shape,Area,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.Properties.VariableNames = ...
        ["Shape","Area (sq km)","Power (dBm)","Power Range (dBm)"];
end

The coveragelevel function returns the coverage level for a latitude-longitude point.

function powerLevels = coveragelevel(lat,lon,GT)

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

    % Find the highest coverage level containing the point
    powerLevels = max(GT.("Power (dBm)")(inContour));

    % Return -inf if the point is not contained within any coverage level
    if isempty(powerLevels)
        powerLevels = -inf;
    end
end

The findNanDelimitedParts helper function finds the indices of the first and last elements of each NaN-separated part of the array x.

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

References

[1] [1] Attachment EngineeringStatement SAT-MOD-20131227-00148. https://fcc.report/IBFS/SAT-MOD-20131227-00148/1031348. Accessed 17 Jan. 2023.