max value inside a circular region in grided data

Hi and thanks in advance.

I have a set of grided wind data (lat, lon, windspeed) that I am processing. Data is grided every 0.11deg between (0-60Deg South and 50-190Deg E) and I want to find the max wind inside a circle of say 30km. How can I do this?

Answer 1

[Voss]

% assuming sample data
lat = -60:0;
lon = 90:160;
wind = randi(10, numel(lat), numel(lon));

% forcing max wind of 100m/s at (lat, lon) = (-25, 115)
[~,latIdx] = min(abs(lat--25));
[~,lonIdx] = min(abs(lon-115));
wind(latIdx, lonIdx) = 100;

% visualize data 
figure
imagesc(lon,lat,wind)
axis xy
colorbar

% define the center of circle
center_lat = -30; 
center_lon = 120; 

R = 6371; % radius of the earth in km

% Haversine formula to calculate distances
dlat = lat - center_lat;
dlon = lon - center_lon;
a = sind(dlat.'/2).^2 + cosd(center_lat) .* cosd(lat.') .* sind(dlon/2).^2;
c = 2 * atan2(sqrt(a), sqrt(1-a));
distances = R * c;

figure
imagesc(lon,lat,distances)
line(lon(lonIdx),lat(latIdx),'Marker','x','Color','r')
axis xy
colorbar

% find points within 3000 km radius
radius_km = 3000;
within_radius = distances <= radius_km;

% get the maximum wind speed within the radius
[max_wind_speed,max_idx] = max(wind(within_radius));
tmp = distances(within_radius);

fprintf('maximum wind speed = %.2f\nat a distance of =  %.2f\n', max_wind_speed, tmp(max_idx));

maximum wind speed = 100.00
at a distance of =  742.97

[Sarvesh]

Thanks @Voss. Just what i needed

[Voss]

You're welcome!

Answer 2

[akshatsood]

Dear @Sarvesh,

I understand that you want to determine the maximum value within a circular region in gridded data. To find the maximum wind speed within a 30km radius for your gridded wind data, you can follow these steps:

Load and Prepare Your Data: Ensure your data is loaded into MATLAB. I am assuming that you have three arrays taking reference from the question: lat, lon and windspeed.

Convert Degrees to Radians: convert your latitude and longitude from degrees to radians using "deg2rad" function.

Calculate Distances: Use the Haversine formula to calculate the distance between your points and the center point of your circle. For more information on Haversone formulae, please refer to the following resource

• https://en.wikipedia.org/wiki/Haversine_formula

Find Maximum Wind Speed: Identify points within the 30km radius and find maximum value amongst these..

% load your data

% define the center of your circle
center_lat = -30; 
center_lon = 100; 

% convert degrees to radians
lat_rad = deg2rad(lat);
lon_rad = deg2rad(lon);
center_lat_rad = deg2rad(center_lat);
center_lon_rad = deg2rad(center_lon);

R = 6371; % radius of the earth in km

% Haversine formula to calculate distances
dlat = lat_rad - center_lat_rad;
dlon = lon_rad - center_lon_rad;
a = sin(dlat/2).^2 + cos(center_lat_rad) .* cos(lat_rad) .* sin(dlon/2).^2;
c = 2 * atan2(sqrt(a), sqrt(1-a));
distances = R * c;

% find points within 30km radius
radius_km = 30;
within_radius = distances <= radius_km;

% get the maximum wind speed within the radius
max_wind_speed = max(windspeed(within_radius));

fprintf('The maximum wind speed within a 30km radius is %.2f\n', max_wind_speed);

I hope this helps.

[Sarvesh]

hi @akshatsood,

Thank you for your response. I tested the code using a set of sample data that I created, and I encountered two issues:

1. Dimension Mismatch: The code stops at 'point a' because the dimensions of the latitude and longitude arrays are not equal (which is the case in my data). How can I address this issue.
2. Incorrect Value Selection: When I adjust the dimensions of the latitude and longitude arrays to match, the final output does not correctly identify the largest value, regardless of the radius size. It seems (and i could be wrong) that the code might be selecting the largest value on the circumference of the circle rather than within it.

% Create sample data
lat = [-60:0];
lon = [90:150]; % error in Haversine form when lon = [90:160]
wind = randi(10, length(lon), length(lat));
wind (35, 30) = 100; % force a sample large value in the data

% Visualize
imagesc(lon, lat, wind)
axis xy
colorbar

% define the center of your circle
center_lat = -25; 
center_lon = 120; 

% convert degrees to radians
lat_rad = deg2rad(lat);
lon_rad = deg2rad(lon);
center_lat_rad = deg2rad(center_lat);
center_lon_rad = deg2rad(center_lon);

R = 6371; % radius of the earth in km

% Haversine formula to calculate distances
dlat = lat_rad - center_lat_rad;
dlon = lon_rad - center_lon_rad;
a = sin(dlat/2).^2 + (cos(center_lat_rad).*cos(lat_rad)).*sin(dlon/2).^2; % point a
c = 2 *atan2(sqrt(a), sqrt(1-a));
distances = R * c;

% find points within 30km radius
radius_km = 3000;
within_radius = distances <= radius_km;

% get the maximum wind speed within the radius
max_wind_speed = max(wind(within_radius));

fprintf('The maximum wind speed within a 30km radius is %.2f\n', max_wind_speed);

[akshatsood]

Dear @Sarvesh

Thank you for response.

To address the issues you are facing, lets break down the solutions:

Dimension Mismatch: The dimension mismatch occurs because the distances are being calculated using latitude and longitude arrays that are not the same size. In such case, code should be modified to calculate distances for each combination of latitude and longitude. This requires creating a grid of latitude and longitude values.

Incorrect Value Selection

To ensure the correct maximum value is selected within the circle, the code needs to compute the distance for each grid point and then check if it is within the specified radius.

Considering the points mentioned above, the code can be modified as follows:

% assuming sample data
lat = -60:0;
lon = 90:150;
wind = randi(10, length(lat), length(lon));

% define the center of circle
center_lat = -25; 
center_lon = 120; 

% convert degrees to radians
lat_rad = deg2rad(lat);
lon_rad = deg2rad(lon);
center_lat_rad = deg2rad(center_lat);
center_lon_rad = deg2rad(center_lon);

R = 6371; % radius of the earth in km

% create a grid of lat and lon values
[lat_grid, lon_grid] = meshgrid(lat_rad, lon_rad);

% Haversine formula to calculate distances
dlat = lat_grid - center_lat_rad;
dlon = lon_grid - center_lon_rad;
a = sin(dlat/2).^2 + cos(center_lat_rad) .* cos(lat_grid) .* sin(dlon/2).^2;
c = 2 * atan2(sqrt(a), sqrt(1-a));
distances = R * c;

% find points within 3000 km radius
radius_km = 3000;
within_radius = distances <= radius_km;

% get the maximum wind speed within the radius
max_wind_speed = max(wind(within_radius));
fprintf('The maximum wind speed within a 3000 km radius is %.2f\n', max_wind_speed);

[Sarvesh]

Hi @akshatsood,

Thank you for response once again.

I understand the first issue but there are still problems with the second one.

I ran the code for lon = 90:150 and the code works fine, i get max_wind_speed = 100 at a minimum radius_km = 727.7 (any lower value than this should not/will not give 100 as output).

When i run the code for lon = 90:160, i should still get the max_winds_speed at the same distance. This does not happen. instead i get max_wind_speed =100 at at minimum radius_km = 3680.76. This should not happen regardless of data gridding.

% assuming sample data
lat = [-60:1:0]';
lon = [90:1:160]';
wind = randi(10, length(lat), length(lon));

% forcing max wind of 100m/s at (lat, lon) = (-25, 115)
latIdx = lat == -25;
lonIdx = lon == 115;
wind(latIdx, lonIdx) = 100;
[val, idx] = max(wind, [],"all"); % check for max wind val

% visualize data 
imagesc(lon, lat, wind)
axis xy
colorbar

% define the center of circle
center_lat = -30; 
center_lon = 120; 

% convert degrees to radians
lat_rad = deg2rad(lat);
lon_rad = deg2rad(lon);
center_lat_rad = deg2rad(center_lat);
center_lon_rad = deg2rad(center_lon);

R = 6371; % radius of the earth in km

% create a grid of lat and lon values
[lat_grid, lon_grid] = meshgrid(lat_rad, lon_rad);
% lat_grid = lat_rad;
% lon_grid = lon_rad;

% Haversine formula to calculate distances
dlat = lat_grid - center_lat_rad;
dlon = lon_grid - center_lon_rad;
a = sin(dlat/2).^2 + cos(center_lat_rad) .* cos(lat_grid) .* sin(dlon/2).^2;
c = 2 * atan2(sqrt(a), sqrt(1-a));
distances = R * c;

% find points within 3000 km radius
radius_km = 3680.76;
within_radius = distances <= radius_km;

% get the maximum wind speed within the radius
max_wind_speed = max(wind(within_radius));
fprintf('maximum wind speed = %.2f\n', max_wind_speed);
fprintf('at a distance of =  %.2f\n', radius_km);

% latlon1=[-25 115];
% latlon2=[-30 120];
% [d1km d2km]=lldistkm(latlon1,latlon2)

[akshatsood]

Dear @Sarvesh

Thank you for the response

I beleive the issue you're experiencing seems to be related to the indexing and calculation of distances within the grid. When you extend the longitude range, the grid size changes, affecting how distances are calculated and which points fall within the specified radius.