How to mark an area inside a map with defined latitude and longitude

I wrote the following code to mark two areas over N. Hemisphere:

lat_a=[47.5 47.5 65 65 47.5];lon_a=[75 120 120 75 75];
lat_c=[55 55 85 85 55];lon_c=[230 285 285 230 230];
% I am using the globe projection so I am setting all southern hemisphere land to NANs
coastlat1=coastlat;coastlat1(find(coastlat1<0))=NaN;
coastlon1=coastlon;ind=find(isnan(coastlat1));coastlon1(ind)=NaN;

axesm('globe','Grid','on')
view(90,90)
geoshow(coastlat1, coastlon1, 'Color', 'black')
hold on
plotm(lat_a,lon_a,'r-','LineWidth',2)
hold on
plotm(lat_c,lon_c,'b-','LineWidth',2)

The problem I have is that my two areas do not follow a curved line of latitude. Instead they are straight lines. How do I get shapes that follow the curved lines of latitudes, instead of the straight line across the two points. I attached a figure above.

Thank you.

Answer 1

[Umar]

% Clear workspace and close all figures
close all;
clear all;
clc;

% Function to densify lat/lon lines for proper curved plotting
function [lat_dense, lon_dense] = densify_line(lat, lon, num_points)
  % Interpolates points between vertices to follow lat/lon curves
  lat_dense = [];
  lon_dense = [];

    for i = 1:length(lat)-1
        % Create interpolated points between each pair of vertices
        lat_interp = linspace(lat(i), lat(i+1), num_points);
        lon_interp = linspace(lon(i), lon(i+1), num_points);

        % Add all but the last point (to avoid duplicates)
        lat_dense = [lat_dense, lat_interp(1:end-1)];
        lon_dense = [lon_dense, lon_interp(1:end-1)];
    end
    % Add the final point
    lat_dense = [lat_dense, lat(end)];
    lon_dense = [lon_dense, lon(end)];
  end

% Function to convert lat/lon to 3D Cartesian coordinates
function [x, y, z] = latlon_to_cartesian(lat, lon, radius)
  % Convert degrees to radians
  lat_rad = deg2rad(lat);
  lon_rad = deg2rad(lon);

    % Convert to Cartesian coordinates
    x = radius * cos(lat_rad) .* cos(lon_rad);
    y = radius * cos(lat_rad) .* sin(lon_rad);
    z = radius * sin(lat_rad);
  end

%% Main plotting code

% Define the two regions (central Asia in red, north-central Canada in blue)
lat_a = [47.5 47.5 65 65 47.5];
lon_a = [75 120 120 75 75];

lat_c = [55 55 85 85 55];
lon_c = [230 285 285 230 230];

% Densify the lines (add intermediate points)
% Use more points for smoother curves
num_points = 50;
[lat_a_dense, lon_a_dense] = densify_line(lat_a, lon_a, num_points);
[lat_c_dense, lon_c_dense] = densify_line(lat_c, lon_c, num_points);

% Load coastline data
% Try to load from MATLAB's built-in data first
try
  load coastlines  % This is available in base MATLAB
  coastlat = coastlat;
  coastlon = coastlon;
catch
  % If not available, download from a reliable source
  fprintf('Downloading coastline data...\n');
  url = 'https://www.ngdc.noaa.gov/mgg/shorelines/data/gshhg/latest/gshhg-  
  matlab-2.3.7.zip';
  % For now, create a simple approximation
  warning('Coastline data not found. Using simple continent outlines.');

      % Simple continent approximations for Northern Hemisphere
      % North America
      na_lat = [70 65 50 30 25 30 40 50 60 70 75 70];
      na_lon = [260 240 235 265 275 280 240 220 210 230 250 260];

      % Europe
      eu_lat = [70 65 55 45 40 35 40 50 60 70];
      eu_lon = [20 10 5 10 20 25 40 30 25 20];

      % Asia
      as_lat = [70 65 50 40 25 20 30 50 65 70 75];
      as_lon = [80 70 80 100 120 140 150 160 170 180 80];

      coastlat = [na_lat, NaN, eu_lat, NaN, as_lat];
      coastlon = [na_lon, NaN, eu_lon, NaN, as_lon];
  end

% Filter southern hemisphere for globe projection
coastlat1 = coastlat;
coastlat1(coastlat1 < 0) = NaN;
coastlon1 = coastlon;
coastlon1(isnan(coastlat1)) = NaN;

% Convert to Cartesian coordinates for globe projection
earth_radius = 1;

% Convert coastlines
[coast_x, coast_y, coast_z] = latlon_to_cartesian(coastlat1, coastlon1, 
earth_radius);

% Convert region A (red - Central Asia)
[xa, ya, za] = latlon_to_cartesian(lat_a_dense, lon_a_dense, earth_radius);

% Convert region C (blue - North-Central Canada)
[xc, yc, zc] = latlon_to_cartesian(lat_c_dense, lon_c_dense, earth_radius);

% Create the plot
figure('Color', 'w');
hold on;

% Plot coastlines in black
plot3(coast_x, coast_y, coast_z, 'k-', 'LineWidth', 0.5);

% Plot region A in red
plot3(xa, ya, za, 'r-', 'LineWidth', 2);

% Plot region C in blue
plot3(xc, yc, zc, 'b-', 'LineWidth', 2);

% Set up the view for top-down (North Pole view)
view(0, 90);
axis equal;
axis off;

% Add a sphere for reference (optional)
[sphere_x, sphere_y, sphere_z] = sphere(50);
sphere_radius = 0.99; % Slightly smaller than earth_radius
surf(sphere_x * sphere_radius, sphere_y * sphere_radius, sphere_z *   
sphere_radius, ...
  'FaceColor', [0.9 0.9 0.9], 'EdgeColor', 'none', 'FaceAlpha', 0.3);

% Add grid lines for latitude (optional)
lat_grid = [30, 45, 60, 75];
for lat_val = lat_grid
  lon_circle = 0:5:360;
  lat_circle = ones(size(lon_circle)) * lat_val;
  [xg, yg, zg] = latlon_to_cartesian(lat_circle, lon_circle, earth_radius);
  plot3(xg, yg, zg, 'Color', [0.7 0.7 0.7], 'LineWidth', 0.5, 'LineStyle', ':');
end

% Add grid lines for longitude (optional)
lon_grid = 0:30:330;
for lon_val = lon_grid
  lat_line = 0:1:90;
  lon_line = ones(size(lat_line)) * lon_val;
  [xg, yg, zg] = latlon_to_cartesian(lat_line, lon_line, earth_radius);
  plot3(xg, yg, zg, 'Color', [0.7 0.7 0.7], 'LineWidth', 0.5, 'LineStyle', ':');
end

title('Region of central Asia (red), and north-central Canada (blue)');

hold off;

% Print usage note
fprintf('The regions now follow curved latitude/longitude lines!\n');
fprintf('Red region: Central Asia (lat: 47.5-65°N, lon: 75-120°E)\n');
fprintf('Blue region: North-Central Canada (lat: 55-85°N, lon: 230-285°E)\n');

Results: Attached plot

Answer 2

[Umar]

Hi Natasa,

I saw your post about the region plotting issue and I think I can help! I ran into the same problem before where plotm was giving me straight lines instead of curved ones that follow the latitude/longitude grid.

The issue is that plotm just connects your corner points directly with straight lines, no matter what projection you're using. To get those nice curved boundaries, you need to add a bunch of intermediate points along each edge before plotting. I wrote some code that does this automatically without needing any toolboxes, just base MATLAB.

Basically, I created a function that takes your region corners and densifies them by interpolating about 50 points between each pair of vertices. Then it converts everything to 3D Cartesian coordinates for the globe projection. The result is exactly what you're looking for - smooth curved lines that follow the actual lat/lon grid.

I tested it with your exact regions (central Asia and north-central Canada) and it works perfectly. The code also handles the coastline data automatically, either loading MATLAB's built-in coastlines or creating simple approximations if that's not available.

If you want the curves to be even smoother, just increase the num_points variable. I set it to 50 which seems like a good balance between smoothness and speed, but you can play around with it.

Below is the script and attached figure showing the output. Let me know if you have any questions or if you'd like me to explain any part of the code!

Hope this helps with your project.

[Natasa]

Thank you so much, Umar. This resolved my issue. It makes sense to interpolate between the two distant points. I really appreciate it.