Customize Buildings for Ray Tracing Analysis

Read Buildings from OpenStreetMap File

Specify the name of an OpenStreetMap file [1] containing data for several city blocks in Chicago.

filename = "chicago.osm";

Read the building parts layer from the file by using the readgeotable function and specifying the Layer name-value argument as "buildingparts". The readgeotable function derives information about the buildings from the file and stores the result in a geospatial table. A geospatial table is a table or timetable with a Shape table variable and attribute table variables. The table represents the building parts using polygon shapes in geographic coordinates.

buildings = readgeotable(filename,Layer="buildingparts");

Add the OpenStreetMap basemap as a custom basemap. View the building footprints over the OpenStreetMap basemap.

basemapName = "osm";
url = "a.tile.openstreetmap.org/${z}/${x}/${y}.png"; 
attribution = "©OpenStreetMap contributors";
addCustomBasemap(basemapName,url,Attribution=attribution)

figure
geobasemap(basemapName)
geoplot(buildings,FaceColor="#808080",FaceAlpha=1)
geotickformat dd

Specify Region of Interest

Specify a ROI. To use a predefined region, specify interactivelySelectROI as false. Alternatively, you can interactively select four points that define a region by specifying interactivelySelectROI as true.

interactivelySelectROI = false;
if interactivelySelectROI
    [latROI,lonROI] = ginput(4);
else
    latROI = [41.8794; 41.8782; 41.8782; 41.8794];
    lonROI = [-87.6293; -87.6293; -87.6277; -87.6277];
end

Create a polygon shape from the ROI. Display the ROI on the same map.

latROI(end+1) = latROI(1);
lonROI(end+1) = lonROI(1);
ROI = geopolyshape(latROI,lonROI);

hold on
geoplot(ROI,FaceColor="m",FaceAlpha=0.2,EdgeColor="m",LineWidth=2)

Extract Buildings Within Region of Interest

Get the latitude and longitude limits of the ROI by querying the bounds of the latitude and longitude coordinates.

[latmin,latmax] = bounds(latROI);
[lonmin,lonmax] = bounds(lonROI);

Get the polygon shapes from the geospatial table. Clip the shapes to the ROI.

shape = buildings.Shape;
clipped = geoclip(shape,[latmin latmax],[lonmin,lonmax]);

Find the shapes that lie entirely or partially within the ROI by using logical indexing. When a polygon shape lies completely outside the ROI, the clipped shape has no coordinate data and the NumRegions property of the shape is 0.

idxInsideROI = clipped.NumRegions > 0;

Create a new geospatial table that contains only the buildings that are entirely or partially within the ROI.

buildingsROI = buildings(idxInsideROI,:);

Display the buildings within the ROI on a new map.

figure
geobasemap(basemapName)
geoplot(buildingsROI,FaceColor="#808080",FaceAlpha=1)

Customize Buildings

Customize three buildings within the ROI by changing the building materials to glass, brick, and metal.

Specify the coordinates of a point within each building that you want to change. To use three predefined points, specify interactivelySelectBuildings as false. Alternatively, you can interactively select three points by specifying interactivelySelectBuildings as true.

interactivelySelectBuildings = false;
if interactivelySelectBuildings
    [latGlass,lonGlass] = ginput(1); 
    [latBrick,lonBrick] = ginput(1); 
    [latMetal,lonMetal] = ginput(1); 
else
    latGlass = 41.8788;
    lonGlass = -87.6288;
    latBrick = 41.8790;
    lonBrick = -87.6283;
    latMetal = 41.8786;
    lonMetal = -87.6283;
end

Create point shapes from the coordinates.

pointGlass = geopointshape(latGlass,lonGlass);
pointBrick = geopointshape(latBrick,lonBrick);
pointMetal = geopointshape(latMetal,lonMetal);

Find the building associated with each point by using a loop. For each building part and each point, determine whether the point is in the building part. If the point is in the building part, update the material and color. Use a dark blue-green color for the glass building, a light red color for the brick building, and a dark gray color for the metal building. If the building part does not contain any of the points, then specify the building material as concrete and the building color as light gray. Specify the colors using hexadecimal color codes.

for row = 1:height(buildingsROI)
    bldg = buildingsROI.Shape(row); 
    if isinterior(bldg,pointGlass)
        buildingsROI.Material(row) = "glass";
        buildingsROI.Color(row) = "#35707E";
    elseif isinterior(bldg,pointMetal)
        buildingsROI.Material(row) = "metal";
        buildingsROI.Color(row) = "#151513";
    elseif isinterior(bldg,pointBrick)
        buildingsROI.Material(row) = "brick";
        buildingsROI.Color(row) = "#AA4A44";
    else 
        buildingsROI.Material(row) = "concrete";
        buildingsROI.Color(row) = "#808080";
    end
end

Import Buildings into Site Viewer

Import the buildings into Site Viewer. Display the buildings over the OpenStreetMap basemap.

viewer = siteviewer(Buildings=buildingsROI,Basemap=basemapName);

View the materials stored in the Site Viewer object by querying the Materials property of the object. By default, ray tracing analysis functions such as raytrace use the materials stored in the MatchedCatalogMaterial table variable.

viewer.Materials

ans=4×2 table
     Material     MatchedCatalogMaterial
    __________    ______________________

    "brick"             "brick"         
    "concrete"          "concrete"      
    "glass"             "glass"         
    "metal"             "metal"         

Perform Ray Tracing Analysis

Create a transmitter site and a receiver site.

tx = txsite(Latitude=41.879090,Longitude=-87.628607);
rx = rxsite(Latitude=41.878410,Longitude=-87.628545);

Create a ray tracing propagation model, which MATLAB® represents using a RayTracing object. Configure the model to find propagation paths with up to three surface reflections. By default, the model uses the shooting and bouncing rays (SBR) method.

pm = propagationModel("raytracing",MaxNumReflections=3);

Calculate the propagation paths and return the result as a cell array of comm.Ray objects. Extract the propagation paths from the cell array. Then, display the transmitter site, the receiver site, and the propagation paths.

rays = raytrace(tx,rx,pm);
rays = rays{1};

show(tx)
show(rx)
plot(rays)

View information about a path by clicking it. The information box displays the interaction materials.

You can also get information about propagation paths by viewing the properties of the comm.Ray objects. View the properties of the ninth path.

ray9 = rays(9)

ray9 = 
  Ray with properties:

      PathSpecification: 'Locations'
       CoordinateSystem: 'Geographic'
    TransmitterLocation: [3×1 double]
       ReceiverLocation: [3×1 double]
            LineOfSight: 0
           Interactions: [1×3 struct]
              Frequency: 1.9000e+09
         PathLossSource: 'Custom'
               PathLoss: 91.7246
             PhaseShift: 0.1253

   Read-only properties:
       PropagationDelay: 3.0369e-07
    PropagationDistance: 91.0454
       AngleOfDeparture: [2×1 double]
         AngleOfArrival: [2×1 double]
        NumInteractions: 3

The NumInteractions property shows that the ninth path has three surface reflections. View the reflection materials for the path by querying the Interactions property of the path.

ray9.Interactions.MaterialName

ans = 
"glass"

ans = 
"metal"

ans = 
"concrete"

[1] The OpenStreetMap file is downloaded from https://www.openstreetmap.org, which provides access to crowd-sourced map data all over the world. The data is licensed under the Open Data Commons Open Database License (ODbL), https://opendatacommons.org/licenses/odbl/.