How to select area precisely based on latitude and longitude?

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.

Answer 1

[KSSV]

REad about inpolygon. If you have the coordinates of the closed polygon you can get the indices lying inside this closed polygon.

[IMC]

I don't have exact coordinates of the polygon. I tried to do it with inpolygon. And in case of inpolygon I selected coordinates just on the base of my assumption (which are obviously not very precise). Is there anyother way to select the coordinates (of the defined area) and then use inpolygon?

[KSSV]

Read about getpts and ginput. These functions can be used to extract points.

[Walter Roberson]

The situation is more complicated. These are geographic regions, so the boundaries are not straight lines and inpolygon() does not apply unless you approximate the curve by sufficiently short line segments.

[IMC]

@Walter RobersonSir can you suggest me how to do this task?

[IMC]

@Image Analyst I didn't get the required output using getpts and ginput. Using ginput and selecting the number of points results in following:

[x,y] = ginput(5)

x =
   -0.0582
    0.0079
    0.0771
   -0.0110
   -0.0645
y =
    0.8322
    0.9234
    0.8731
    0.7724
    0.8322

What I want is lat lon values of ploygon rather than above values.

[KSSV]

How you have plotted black lines in the plot?

If you have x, y points use inpolygon.

Answer 2

[Image Analyst]

Not sure what you mean. Is either of these what you want to do?

1. How about drawpolygon() for interactively selecting/drawing the polygon region?
2. 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.

[IMC]

In my figure every pixel has a latitude and longitude value. I just want to study area specified in polygon. Secondly, I don't have the exact lat lon values of this polygon (and I need to be very precise about it). My question is what should I do to get the lat lon values of this polygon? Later I'll just use this (polygon lat lon values) for indexing purpose. (For example, to study cloud temperatue (CT) within this polygon area i.e. CT(lat)).

[Image Analyst]

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).

[IMC]

Lat and Lon limit I set for the figure:

h1=axesm('mercator','MapLatLimit',[15 55],'MapLonLimit',[80 150]);

Secondly, in this hdf image which I have displayed every pixel has its lat lon value. I just read the parameter (Cloud temperature(CT)) and lat, lon using hdfread. And then made the figure.

geoshow(h1,lat,lon,CT,'DisplayType','texturemap')

Min Lat of the image: 33.1063

Max Lat of the image: 54.5416

Min Lon of the image: 95.9691

Max Lon of the image: 130.7437

Figure is attached below. Thank you

Answer 3

[Image Analyst]

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

[Walter Roberson]

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)

[IMC]

@Image Analyst Thank you! I'll try this.