max value inside a circular region in grided data

3 views (last 30 days)
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?

Accepted Answer

Voss
Voss on 16 Aug 2024
% 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

More Answers (1)

akshatsood
akshatsood on 15 Aug 2024
Edited: akshatsood on 15 Aug 2024
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
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.
  4 Comments
Sarvesh
Sarvesh on 16 Aug 2024
Edited: Sarvesh on 16 Aug 2024
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
akshatsood on 16 Aug 2024
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.

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!