Where are you flying over?

저자 Jiro Doke, May 20, 2016

1 회 조회 (최근 30일) | 0 좋아요 | 0 댓글

Jiro's pick this week is CityName by Richard Moore.

When I used to travel to various places in the U.S. for work, I would wonder which cities I was flying over during the flight. On some planes, they have seat back monitors that show the flight path with names of cities the plane is flying over.

With Richard's CityName, I can do the same. I simply call the function like this

[city, country, dist] = CityName(lat, lon)

and it returns the closest city to the latitude and logitude location based on the GeoNames 5000 cities. It also returns the distance (in km) from the city center.

Lets say that I am flying from Boston, Massachusetts to Los Angeles, California.

boston = [42.36 -71.06];
los = [34.052 -118.244];

We can plot the map of the U.S. using the Mapping Toolbox.

ax = usamap('conus');
states = shaperead('usastatehi', 'UseGeoCoords', true,...
  'Selector', {@(name) ~any(strcmp(name,{'Alaska','Hawaii'})), 'Name'});
geoshow(ax, states, 'DisplayType', 'polygon', 'FaceColor', [0.5 1 0.5])
framem off, gridm off, mlabel off, plabel off

Next, we can calculate the great circle track between the two cities, again using track2 from the Mapping Toolbox. Yes, this is not the path a plane would take, but for the purpose of this post, this will do.

[lat,lon] = track2(boston(1),boston(2),los(1),los(2));

Now, we can write a loop for each coordinate to find the closest city within 10 km of the path.

% Figure out view limits for animation purpose
sLim = [1386896     3288804     4332357     5404342];
eLim = [-2438739     -536831     3335946     4407931];
limits = [linspace(sLim(1),eLim(1))',linspace(sLim(2),eLim(2))',...
    linspace(sLim(3),eLim(3))',linspace(sLim(4),eLim(4))'];

right = true;   % Used for text placement
for id = 1:100
    % Plot a segment for the path
    if id < 100
        plotm(lat(id:id+1),lon(id:id+1),'Color','r','LineWidth',2)
    end

    % Set axis limits
    axis(limits(id,:))

    % Find closest city
    [c,~,dist] = CityName(lat(id),lon(id));

    % If the distance is within 10 km, display name of the city
    if dist < 10
        c = strtrim(c);

        % Place point
        plotm(lat(id),lon(id),'Color','b','Marker','o','MarkerFaceColor','b')

        % Display name of city. Alternate placing text on right and left.
        if right
            textm(lat(id),lon(id),c,...
                'HorizontalAlignment','left','VerticalAlignment','top',...
                'FontSize',12,'FontWeight','bold','Color','r')
        else
            textm(lat(id),lon(id),c,...
                'HorizontalAlignment','right','VerticalAlignment','bottom',...
                'FontSize',12,'FontWeight','bold','Color','r')
        end
        right = ~right;
    end
    drawnow
end