Как определить круговые границы в MATLAB

Question

Я работаю с кодом MATLAB, и я хотел бы определить граничные условия.
У меня есть данные широты молнии, хранящиеся в столбце 3, и данные о долготе, хранящиеся в столбце 4.
Кроме того, у меня есть координата глаза тропического циклона.
Теперь у меня есть 3 области, которые мне нужно отделить, а именно:

1. eyewall (0 - 60 км)
2. внутренняя rainband (60-180 км) и,
3. наружный rainband (180 - 500 км). от глаз циклона.

Я преобразовал расстояния в терминах степеней широты, чтобы облегчить дело.
То, что я хотел бы сделать, - подсчитать случаи возникновения молний отдельно в этих трех областях циклона, когда глаз находится на заданной широте и долготе. Время хранится во втором столбце.
Я попытался использовать математический эквивалент окружности с центром (а, б), соответствующий центром циклона и г как расстояние от центра:
(Xa) + (Yb) = г , с й, представляющим столбцом долготы 4 и у, представляющего столбец широты 3.
я пришел с последующим кодом для конкретного случая, но это не похоже на работу или скомпилировать:

M = A20121228; 
lat = -10.6; 
lon = 161.5; 
hr = 0; 

if ((M(:, 4) - lon)^2 + (M(:, 3) - lat)^2) <= 0.541^2 
row_idx1 = (hr == M(:, 2) & lat == M(:, 3) & lon == M(:, 4)); 
filtered_M1 = M(row_idx1, :); 
eyewall = filtered_M1; 
end  

if ((M(:, 4) - lon)^2 + (M(:, 3) - lat)^2) > 0.541^2 && ((M(:, 4) - lon)^2 + (M(:, 3) - lat)^2) <= 1.622^2 
row_idx2 = (hr == M(:, 2) & lat == M(:, 3) & lon == M(:, 4)); 
filtered_M2 = M(row_idx2, :); 
inner = filtered_M2; 
end 

if ((M(:, 4) - lon)^2 + (M(:, 3) - lat)^2) > 1.622^2 && ((M(:, 4) - lon)^2 + (M(:, 3) - lat)^2) <= 4.505^2 
row_idx3 = (hr == M(:, 2) & lat == M(:, 3) & lon == M(:, 4)); 
filtered_M3 = M(row_idx3, :); 
outer = filtered_M3; 
end 

Answer 1

Вы должны сформировать свои геодезические данные в терминах «километр» вместо конвертируя диапазоны (0, 60, 180, 500) с точки зрения степеней, потому что работа на «километр» устройства более удобна в вашем случае. Вы можете использовать функцию geodetic2ned(), чтобы преобразовать ваши геодезические данные в локальные декартовы NED (Nort - East - Down -> они являются осью Navigation Frame). При использовании geodetic2ned(), если вы введете свои lat = -10.6 и lon = 161.5 (я назвал их lat0 и lon0) в качестве аргумента для ссылки, они станут источником навигационной рамки. Поэтому N, E (являются выходами функции geodetic2ned(), проверьте код ниже) будут оценены по вашему адресу lat и lon (иными словами lat0 и lon0). К таким образом, вы можете написать непосредственно sqrt(N.^2 + E.^2) вместо sqrt(M(:, 4) - lon)^2 + (M(:, 3) - lat)^2) (кстати, вы должны использовать Element-Wise Power)

Теперь проверьте код ниже:

SPHEROID = referenceEllipsoid('wgs84', 'km'); % // SPHEROID Reference 
hr   = 0; 
lat0  = -10.6; 
lon0  = 161.5; 
[N, E, D] = geodetic2ned(M(:, 3), M(:, 4), -M(:, 2), lat0, lon0, hr, SPHEROID); 
range  = sqrt(N.^2 + E.^2); 
indEyewall = find(hr == D & range >= 0 & range < 60); 
indInner = find(hr == D & range >= 60 & range < 180); 
indOuter = find(hr == D & range >= 180 & range < 500); 
eyewall = M(indEyewall, :); 
inner  = M(indInner, :); 
outer  = M(indOuter, :); 

Дело, я включаю hr == D часть для find() функции, потому что вы пишете в вашем коде, но я не уверен, если вы включите его. Вы можете получить логическое значение false от find() функций, поэтому, если бы я был вами, я бы позаботился о его включении.

Надеюсь, это было бы полезно. Удачи ..