MATLAB: How to access data in structure by finding one of its elements

cellfindsearchstructstructure

I am making project using data from map. I have 2 structures. "ways" and "nodes". Nodes contain ID, Latitude and Longitude. Ways contain information about "which nodes (IDs) represent this road (they are points marking street)" and data I want to retrieve. Structures are following:

 nodes = 
      version: 0.6000
    generator: 'Overpass API'
        osm3s: [1x1 struct]
     elements: {37952x1 cell}

In elements there are cells containing 1×1 struct like that:

 nodes.elements{1,1} = 
    type: 'node'
      id: 30788202
     lat: 54.4719
     lon: 18.4729

———————

Ways are similar:

 ways = 
      version: 0.6000
    generator: 'Overpass API'
        osm3s: [1x1 struct]
     elements: [247x1 struct]

Inside:

 ways.elements = 
 247x1 struct array with fields:
    type
    id
    nodes
    tags

For example:

 ways.elements(1) = 
     type: 'way'
       id: 4941052
    nodes: [6x1 double]
     tags: [1x1 struct]

—————————

Now. I have GPS positions from my ride (a lot of them) – lat and lon. I create small box around each of my positions and then search all the nodes to check if any of them is inside my box. If I find that I am near the node, I search all the ways to find what road is it (I check ways.elements(i).nodes if they contain the node near me). When found I take necessary data from tags.

————————-

Code looks like that:

 diff = 0.1e-4; %to create small box around my position
 %initialization
 speedlimit = zeros(length(GPSinterpol),1);
 previousNode = 0;                               %just to make sure they are different

 nodeID = -1;                                    %just to make sure they are different
 for k=1:length(GPSinterpol)                     %for all positions I have been at
    position = GPSinterpol(k,:);
    positionplus = position + diff;              %create upper right corner of the box
    positionminus = position - diff;             %create lower left corner of the box
    for i = 1:length(nodes.elements)             %search all nodes
                                                 %check if any node is inside my box
        if all([nodes.elements{i,1}.lat nodes.elements{i,1}.lon] > positionminus) && all([nodes.elements{i,1}.lat nodes.elements{i,1}.lon] < positionplus)
            nodeID = nodes.elements{i,1}.id;
            if nodeID ~= previousNode
                %find way corresponding our node
                for j = 1:length(ways.elements)                                           %search all ways
                    if any(ways.elements(j,1).nodes == nodeID)                            %way found
                        speedlimit(k) = str2double(ways.elements(j,1).tags.maxspeed);     %get speedlimit from tags
                        break
                    end
                end
            end
            break
        end
    end
    if nodeID == previousNode
        speedlimit(k) = speedlimit(k-1);
    end
    previousNode = nodeID;
 end

My question is how to do this better ;). Using loops is extremely slow for that amount of data and I want to use it on even bigger. I hope there are better ways to search structures

Best Answer

First, don't use diff as a variable name since it will stop you using the very useful diff function.

Secondly, it's always dangerous to use length on 2D data, since it returns the length of the largest dimension which you assume is the number of rows. Maybe one day, you'll have only one GPS position, and length will return the number of columns causing failure later on the code. size(GPSinterpol, 1) is guaranteed to return the number of rows no matter what, so why not use it.

Thirdly, from your code, it looks like a node can only belong to one way, or if a node belongs to more than one way, you always select the last way that include that node. So, instead of doing the lookup of which way a node belongs to in a loop, you could do it just once for all the nodes at the start. Personally, I would add the way id to the node structure, but you could also use a containers.map or a simple vector;

waynodes = {way.elements.nodes}  %concatenate all the way nodes into a cell array
wayids = repelem([way.elements.id], cellfun(@numel, waynodes));  %replicate the way ids to match the nodes
waynodes = cell2mat(waynodes);   %convert cell array into column vector
%store as plain matrix that can be looked up:
nodewaymap = [waynodes, wayids'];
%looking up the way id for a node is
nodewaymap(nodewaymap(:, 1) == nodeID), 2)
%store as containers.Map
nodewaymap = containers.Map(waynodes, wayids);  %may fail if node belongs to more than one way
%looking up the way id for a node is simply
nodewaymap(nodeID)
%or add as new field of node structure
nodesstruct = [nodes.elements{:}];  %the cell array complicates everything
nodesid = [nodesstruct.id];
[~, order] = ismember(nodesid, waynodes);  %I'm assuming all nodes belong to a way
wayids = num2cell(wayids(order));
[nodesstruct.wayid] = deal(wayids{:});  %add new field in the right order
%looking up is then straightforward:
nodesstruct.wayid

As you can see in the comment above, the storage of the nodes as scalar structure in a cell array, complicates things. You'll be better off converting that cell array into a an array structure:

nodesstruct = [nodes.elements{:}];

Finally, you can indeed get rid of the node searching loop:

positiondiff = 0.1e-4
positionplus = GPSinterpol + positiondiff;
positionminus = GPSinterpol + positionminus;
nodeslatlon = [cat(3, nodesstruct.lan), cat(3, nodesstruct.lon)];
%in R2016b or later:
ismatch = squeeze(all(nodeslatlon > positionminus & nodeslatlon < positionplus));
%versions < R2016b:
ismatch = squeeze(all(bsxfun(@gt, nodeslatlon, positionminus) & bsxfun(@lt, nodeslatlon, positionplus)));
%rows of ismatch are the GPSlocation, columns are the nodes. A 1 indicates that location and node match, so
[~, nodeindex] = max(ismatch, 2);  %returns the index of the node that correspond to locations

And depending on the method you used to store the wayid:

%as matrix:
[~, rowindex] = ismember(nodeindex, nodewaymap(:, 1));
gspwayid = nodewaymap(rowindex, 2);
%as map
[~, mapindex] = ismember(nodeindex, nodewaymap.keys);
gpswayid = [nodewaymap.values{mapindex}];
%as field
gpswayid = [nodesstruct(nodeindex).wayid];

And converting to speed:

[~, wayindex] = ismember(gpswayid, [way.elements.id]);
tags = [way.elements(wayindex).tags];
speedlimit = str2double({tags.maxspeed});

And the above just made me realise that storing the structure index of the ways instead of the id, would avoid that last code to ismember.

Related Solutions

MATLAB: How to generate a mesh with quadratic tetrahedral elements using function ‘generateMesh’ whose boundary nodes follow the underlying curved geometry

This feature is unavailable as per MATLAB R2020b, as the high order nature of the elements is not based on the fact that the boundary faces of these elements match the geometry itself, but that they rather offer more nodes for a quadratic order approximation of the unkown field, which is standard in the industry.

Herein it is presented a simplified demo based on a tetrahedral discretization of a sphere. The workflow provided below allows for modifying the coordinates of these additional nodes that are located on the corresponding boundary face of the tetrahedral elements, such that all the underlying boundary nodes lie exactly on the spherical boundary shape. Moreover, this demo only shows how the corresponding nodal coordinates of the of the first element in the mesh that happens to be on the boundary of the sphere can be accordingly modified, such that its underlying nodes lie on the sperical boundary. This demo can be extended with a loop over all the elements, whereby the corresponding boundary elements also need to be identified, before their boundary nodes can be adjusted and it is only meant to serve as a demonstrator.

%% Create a sphere and a high order mesh
Radius = 5;
sp = multisphere(Radius);
m1 = createpde;
m1.Geometry = sp;
Hmax = 2;
mesh = generateMesh(m1,'Hmax',Hmax);
pdeplot3D(m1);
% Plot the original mesh that results when using function 'generateMesh'
% hold on
% for i = 1:length(m1.Mesh.Elements(:, 1))
%     nodeID = m1.Mesh.Elements(i, 1);
%     plot3(m1.Mesh.Nodes(1, nodeID), m1.Mesh.Nodes(2, nodeID), m1.Mesh.Nodes(3,nodeID), 'ro');
% end
%% Modify the nodal locations according to the spherical shape
% define the center of the sphere
center = [0; 0; 0];
% get the elements and the nodes of the mesh
els = m1.Mesh.Elements;
nds = m1.Mesh.Nodes;
% initialize auxiliary variables
displVct = zeros(3, length(els(:,1)));
displ = zeros(length(els(:,1)), 1);
% loop over all the nodes of the first element
for i = 1:length(els(:,1))
  % get the nodal ID
  nodeID = els(i,1);
  % compute the displacement vector and its magnitude from the center of
  % the sphere
  distCenterVct = nds(:,nodeID);
  distCenter = norm(distCenterVct);
  % compute the spherical coordinates of the node
  [azimuth, elevation, ~] = ...
    cart2sph(nds(1,nodeID), nds(2,nodeID), nds(3,nodeID));
  % compute the position vector of the point on the exact sphere
  vec = [Radius*cos(elevation)*cos(azimuth)
    Radius*cos(elevation)*sin(azimuth)
    Radius*sin(elevation)];
  % compute the necessary displacement to bring that node to the exact
  % spherical shape
  displVct(:, i) = vec - distCenterVct;
  displ(i, 1) = norm(displVct(:, i));
end
% sort the nodes by their displacement to the exact sphere
[displSorted, idx] = sort(displ);
% The three first idx(1:3) closest nodes are the vertex nodes of the 
% tetrahedral that lie on the sphere (due to approximation errors the 
% corresponding displacements won't be exactly zero)
%

% The three ones afterwards idx(4:6) are the high order nodes that lie
% within the edges of the tetrahedral at the boundary face of the 
% tetrahedral
%
% The rest of the nodes idx(7:10) are interior to the tetrahedral and have
% no influence on the outer shape
nodesToBeMoved = idx(4:6);
% loop over all nodes whose nodal locations are to be updated
for i = 1:length(nodesToBeMoved)
  % get the nodal IDs
  nodeID = els(nodesToBeMoved(i, 1),1);
  % update the nodal locations by the necessary displacements to achieve
  % that they lie on the sphere
  nds(:, nodeID) = nds(:, nodeID) + displVct(:, nodesToBeMoved(i, 1));
end
%% create a new pde problem and assign the new elements and nodes with modified locations to it
m2 = createpde;
geometryFromMesh(m2,nds,els);
figure
pdeplot3D(m2);
hold on
for i = 1:length(m2.Mesh.Elements(:, 1))
    nodeID = m2.Mesh.Elements(i, 1);
    plot3(m2.Mesh.Nodes(1, nodeID), m2.Mesh.Nodes(2, nodeID), m2.Mesh.Nodes(3,nodeID), 'ro');
end

Running the aforementioned code snippet results in having the boundary nodes of the first element, that also happens to lie on the boundary, lying on the exact spherical shape, please see the following screenshot,

Please note that the tetrahedrals are herein only visualized by their 4 corner nodes and not by the additional ones that are added by the quadratic approximation.

MATLAB: How to convert a cell array of structures to array of structures

As long as the structures in the cells have all the same fields and that these structures in the cells are either scalar or arrays with the same number of rows:

nodedata = [nodeinfo.NodeData{:}];
%or
nodedata = cell2mat(nodeinfo.NodeData);

Best Answer

Related Solutions

MATLAB: How to generate a mesh with quadratic tetrahedral elements using function ‘generateMesh’ whose boundary nodes follow the underlying curved geometry

MATLAB: How to convert a cell array of structures to array of structures

Related Question