How to select area precisely based on latitude and longitude?
22 views (last 30 days)
Show older comments
Hello everyone,
I need to do analysis for just this part (black in figure 1). I want to be very precise regarding area selection. Is there anyway I can select area precisely based on latitude and longitude? The data I plotted is from hdf file. Thank you.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/674948/image.jpeg)
0 Comments
Answers (3)
Image Analyst
on 5 Jul 2021
Not sure what you mean. Is either of these what you want to do?
- How about drawpolygon() for interactively selecting/drawing the polygon region?
- Or do you mean that you have a list of N (lat, lon) coordinates already somehow and you want to outline or colorize inside that list on your displayed image.
3 Comments
Image Analyst
on 5 Jul 2021
Somehow, some way, you need to specify the latitude and longitude. Since you don't seem to have a list of them in advance, and you don't seem to want to draw them with drawpolygon(), then how to you want to specify it?
If you have an image (not a figure or screenshot) of the mostly white map with some parts colorized in a polygon, then we can detect which pixels are brightly colored and then get the outline of that. Is that what you want to do? If so, please upload the actual digital image you displayed (not the screenshot with axes, tick marks and labels, color bar, etc.) -- I just want the bare image alone. Also I need the variable that contains the list or formula for giving the lat and lon of every pixel. If you don't have it at least tell us the lat of line 1 (? N) and the last line (15N), and the longitude of the left column (80 maybe???) and the right column (150E I'd guess).
Image Analyst
on 5 Jul 2021
Try this. It asks you for the limits of lat and lon, then spatially calibrates the image with a function that turns (row, column) or (y, x) into latitidue and longitude. It then finds the colorized polygon region and labels the limits of the bounding box with the lat, lon. Then it prints out the lat and lon for every pixel in the polygon region to the command window.
% Demo by Image Analyst, June, 2021.
clc; % Clear the command window.
close all; % Close all figures (except those of imtool.)
clearvars;
workspace; % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 16;
fprintf('Beginning to run %s.m ...\n', mfilename);
%-----------------------------------------------------------------------------------------------------------------------------------
% Read in image.
folder = [];
baseFileName = 'map.jpeg';
fullFileName = fullfile(folder, baseFileName);
% Check if file exists.
if ~isfile(fullFileName)
% The file doesn't exist -- didn't find it there in that folder.
% Check the entire search path (other folders) for the file by stripping off the folder.
fullFileNameOnSearchPath = baseFileName; % No path this time.
if ~exist(fullFileNameOnSearchPath, 'file')
% Still didn't find it. Alert user.
errorMessage = sprintf('Error: %s does not exist in the search path folders.', fullFileName);
uiwait(warndlg(errorMessage));
return;
end
fullFileName = fullFileNameOnSearchPath;
end
rgbImage = imread(fullFileName);
[rows, columns, numberOfColorChannels] = size(rgbImage)
% Display the test image full size.
subplot(2, 2, 1);
imshow(rgbImage, []);
axis('on', 'image');
caption = sprintf('Original Image : "%s"', baseFileName);
title(caption, 'FontSize', fontSize, 'Interpreter', 'None');
drawnow;
hp = impixelinfo(); % Set up status line to see values when you mouse over the image.
% Set up figure properties:
% Enlarge figure to full screen.
hFig1 = gcf;
hFig1.Units = 'Normalized';
hFig1.WindowState = 'maximized';
% Get rid of tool bar and pulldown menus that are along top of figure.
% set(gcf, 'Toolbar', 'none', 'Menu', 'none');
% Give a name to the title bar.
hFig1.Name = 'Demo by Image Analyst';
%--------------------------------------------------------------------------------------------------------
% Find the axes box.
axesBoxMask = rgbImage(:, :, 1) < 128;
% Fill holes
axesBoxMask = imfill(axesBoxMask, 'holes');
% Take largest blob only
axesBoxMask = bwareafilt(axesBoxMask, 1);
subplot(2, 2, 2);
imshow(axesBoxMask);
grid on;
title('Axes Mask', 'FontSize', fontSize);
% Get the bounding box
[r, c] = find(axesBoxMask);
row1 = min(r);
row2 = max(r);
col1 = min(c);
col2 = max(c);
%================================================================================================================================
% Get latitude and longitude coordinates.
% Ask user for two floating point numbers.
defaultValue = {'60', '15'};
titleBar = 'Enter latitude values';
userPrompt = {'Enter top latitude : ', 'Enter bottom latitude: '};
caUserInput = inputdlg(userPrompt, titleBar, 1, defaultValue);
if isempty(caUserInput),return,end % Bail out if they clicked Cancel.
% Convert to floating point from string.
lat1 = str2double(caUserInput{1})
lat2 = str2double(caUserInput{2})
% Check usersValue1 for validity.
if isnan(lat1)
% They didn't enter a number.
% They clicked Cancel, or entered a character, symbols, or something else not allowed.
% Convert the default from a string and stick that into usersValue1.
lat1 = str2double(defaultValue{1});
message = sprintf('I said it had to be a number.\nTry replacing the user.\nI will use %.2f and continue.', lat1);
uiwait(warndlg(message));
end
% Do the same for usersValue2
% Check usersValue2 for validity.
if isnan(lat2)
% They didn't enter a number.
% They clicked Cancel, or entered a character, symbols, or something else not allowed.
% Convert the default from a string and stick that into usersValue2.
lat2 = str2double(defaultValue{2});
message = sprintf('I said it had to be a number.\nTry replacing the user.\nI will use %.2f and continue.', lat2);
uiwait(warndlg(message));
end
%================================================================================================================================
% Get longitude coordinates.
% Ask user for two floating point numbers.
defaultValue = {'80', '150'};
titleBar = 'Enter longitude values';
userPrompt = {'Enter left longitude : ', 'Enter right longitude: '};
caUserInput = inputdlg(userPrompt, titleBar, 1, defaultValue);
if isempty(caUserInput),return,end % Bail out if they clicked Cancel.
% Convert to floating point from string.
lon1 = str2double(caUserInput{1})
lon2 = str2double(caUserInput{2})
% Check usersValue1 for validity.
if isnan(lon1)
% They didn't enter a number.
% They clicked Cancel, or entered a character, symbols, or something else not allowed.
% Convert the default from a string and stick that into usersValue1.
lon1 = str2double(defaultValue{1});
message = sprintf('I said it had to be a number.\nTry replacing the user.\nI will use %.2f and continue.', lon1);
uiwait(warndlg(message));
end
% Do the same for usersValue2
% Check usersValue2 for validity.
if isnan(lon2)
% They didn't enter a number.
% They clicked Cancel, or entered a character, symbols, or something else not allowed.
% Convert the default from a string and stick that into usersValue2.
lon2 = str2double(defaultValue{2});
message = sprintf('I said it had to be a number.\nTry replacing the user.\nI will use %.2f and continue.', lon2);
uiwait(warndlg(message));
end
%================================================================================================================================
% Get colorized region.
[regionMask, maskedRGBImage] = createMask(rgbImage);
% Fill holes
regionMask = imfill(regionMask, 'holes');
% Take largest 2 blob only - the colorized region and the color bar (in case the color bar is larger).
regionMask = bwareafilt(regionMask, 2);
% Take the left most one. It will have label 1
labeledImage = bwlabel(regionMask);
regionMask = ismember(labeledImage, 1); % Extract the region labeled #1.
% Display the image.
subplot(2, 2, 3);
imshow(regionMask, []);
axis('on', 'image');
title('Region Mask', 'FontSize', fontSize, 'Interpreter', 'None');
drawnow;
hp = impixelinfo(); % Set up status line to see values when you mouse over the image.
% Get a formula for converting column into latitude and longitude
% lat = a * row + b
latCoefficients = polyfit([row1, row2], [lat1, lat2], 1)
lonCoefficients = polyfit([col1, col2], [lon1, lon2], 1)
% Get coordinates of the polygon regionMask
[rp, cp] = find(axesBoxMask);
% Get the bounding box
prow1 = min(rp);
prow2 = max(rp);
pcol1 = min(cp);
pcol2 = max(cp);
markerSize = 25;
hold on;
% Label upper left corner.
[lat, lon] = GetLatAndLon(prow1, pcol1, latCoefficients, lonCoefficients)
str = sprintf(' (lat, lon) =\n (%.1f, %.1f)', lat, lon);
plot(pcol1, prow1, 'r+', 'MarkerSize', markerSize, 'LineWidth', 2);
text(pcol1, prow1, str, 'Color', 'r', 'FontSize', 10, 'FontWeight', 'bold');
% Label lower left corner.
[lat, lon] = GetLatAndLon(prow2, pcol1, latCoefficients, lonCoefficients)
str = sprintf(' (lat, lon) =\n (%.1f, %.1f)', lat, lon);
plot(pcol1, prow2, 'r+', 'MarkerSize', markerSize, 'LineWidth', 2);
text(pcol1, prow2, str, 'Color', 'r', 'FontSize', 10, 'FontWeight', 'bold');
% Label upper right corner.
[lat, lon] = GetLatAndLon(prow1, pcol2, latCoefficients, lonCoefficients)
str = sprintf(' (lat, lon) =\n (%.1f, %.1f)', lat, lon);
plot(pcol2, prow1, 'r+', 'MarkerSize', markerSize, 'LineWidth', 2);
text(pcol2, prow1, str, 'Color', 'r', 'FontSize', 10, 'FontWeight', 'bold');
% Label lower right corner.
[lat, lon] = GetLatAndLon(prow2, pcol2, latCoefficients, lonCoefficients)
str = sprintf(' (lat, lon) =\n (%.1f, %.1f)', lat, lon);
plot(pcol2, prow2, 'r+', 'MarkerSize', markerSize, 'LineWidth', 2);
text(pcol2, prow2, str, 'Color', 'r', 'FontSize', 10, 'FontWeight', 'bold');
drawnow;
% Print out the latitude and longitude for every point in the region mask.
for k = 1 : length(rp)
x = cp(k);
y = rp(k);
[lat, lon] = GetLatAndLon(y, x, latCoefficients, lonCoefficients);
fprintf('Point (row, col) = (%d, %d) is at (lat, lon) = (%f, %f).\n', ...
y, x, lat, lon);
end
fprintf('Done running %s.m\n', mfilename);
msgbox('Done. Check command window for coordinates.');
%==============================================================================================
function [BW,maskedRGBImage] = createMask(RGB)
%createMask Threshold RGB image using auto-generated code from colorThresholder app.
% [BW,MASKEDRGBIMAGE] = createMask(RGB) thresholds image RGB using
% auto-generated code from the colorThresholder app. The colorspace and
% range for each channel of the colorspace were set within the app. The
% segmentation mask is returned in BW, and a composite of the mask and
% original RGB images is returned in maskedRGBImage.
% Auto-generated by colorThresholder app on 05-Jul-2021
%------------------------------------------------------
% Convert RGB image to chosen color space
I = rgb2hsv(RGB);
% Define thresholds for channel 1 based on histogram settings
channel1Min = 0.000;
channel1Max = 0.939;
% Define thresholds for channel 2 based on histogram settings
channel2Min = 0.328;
channel2Max = 1.000;
% Define thresholds for channel 3 based on histogram settings
channel3Min = 0.000;
channel3Max = 1.000;
% Create mask based on chosen histogram thresholds
sliderBW = (I(:,:,1) >= channel1Min ) & (I(:,:,1) <= channel1Max) & ...
(I(:,:,2) >= channel2Min ) & (I(:,:,2) <= channel2Max) & ...
(I(:,:,3) >= channel3Min ) & (I(:,:,3) <= channel3Max);
BW = sliderBW;
% Initialize output masked image based on input image.
maskedRGBImage = RGB;
% Set background pixels where BW is false to zero.
maskedRGBImage(repmat(~BW,[1 1 3])) = 0;
end
% Return the longitude and latitude given a row and column (y and x).
function [lat, lon] = GetLatAndLon(prow1, pcol1, latCoefficients, lonCoefficients)
x = pcol1;
y = prow1;
lon = polyval(lonCoefficients, x);
lat = polyval(latCoefficients, y);
end
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/675253/image.png)
2 Comments
Walter Roberson
on 5 Jul 2021
This is not directly usable for the situation, unfortunately.
NetCDF and HDF files often contain nonlinear grids, such as having been measured at equal physical distances but coordinates given in latitude and longitude .
The upper edge of the user's latitude range is 2149 km long (Russia, north of Mongolia) .
The lower edge of the user's latitude range is 3128 km long (China to South Korea)
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!