Calculate pressure gradient between two areas?

4 visualizzazioni (ultimi 30 giorni)
Shayma Al Ali
Shayma Al Ali il 29 Nov 2020
Risposto: Hassaan il 3 Gen 2024
So I have two variables, zonal1 and zonal2. Each variable contains a list of latitude, longitude, and pressures. I need to calculate the pressure gradient which is . The dimensions for zonal1 is 1x85 while zonal2 is 1x168. My question is how do I calculate the distance between the two variables? Do I create a 3D meshgrid for both variables? I'm just kind of confused. My code is listed below:
%Cd to Zonal 1 area
cd /Users/shaymaalali/Desktop/Physical_Oceanography/cchdo_320620170820
%Read Files
direc=dir('*.nc');
zonal1=[];
for i=1:length(direc)
zonal1.lat=ncread(direc(i).name,'latitude'); %get latitude
zonal1.lon=ncread(direc(i).name,'longitude'); %get longitude data
zonal1.time=ncread(direc(i).name,'time'); %get time
zonal1.pressure=ncread(direc(i).name,'pressure'); %get pressure
zonal1.temp=ncread(direc(i).name,'temperature'); %get temp
zonal1.salinity=ncread(direc(i).name,'salinity'); %get salinity
end
%Calculate density
zonal1.density=[];
for j=1:length(zonal1)
%cd pathway
[zonal1.density,alpha, beta]=gsw_rho_alpha_beta_CT_exact(zonal1.salinity,zonal1.temp,zonal1.pressure); %get density from TEOS-10 function
end
%Calculate depth from Pressure
zonal1.depth=[];
for i=1:length(zonal1)
zonal1.depth=Depth(zonal1.pressure,zonal1.lat);
end
%%
%Cd to Zonal 2 area
cd /Users/shaymaalali/Desktop/Physical_Oceanography/cchdo_74EQ20151206
%ZipFiles
unzip('1_a05_74EQ20151206_nc_ctd.zip');
unzip('0_a05_74EQ20151206_ct1.zip');
direc=dir('*.nc');
zonal2=[]; %create variable
for i=1:length(direc)
zonal2.lat=ncread(direc(i).name,'latitude'); %get latitude
zonal2.lon=ncread(direc(i).name,'longitude'); %get longitude data
zonal2.time=ncread(direc(i).name,'time'); %get time
zonal2.pressure=ncread(direc(i).name,'pressure'); %get pressure
zonal2.temp=ncread(direc(i).name,'temperature'); %get temp
zonal2.salinity=ncread(direc(i).name,'salinity'); %get salinity
end
zonal2.density=[];
%pathway='/Users/shaymaalali/Desktop/Physical_Oceanography/gsw_matlab_v3_06_12';
for j=1:length(zonal2)
%cd pathway
[zonal2.density,alpha, beta]=gsw_rho_alpha_beta_CT_exact(zonal2.salinity,zonal2.temp,zonal2.pressure); %get density from TEOS-10 function
end
%calculate depth
zonal2.depth=[];
%Calculate depth from pressure
for i=1:length(zonal2)
zonal2.depth=Depth(zonal2.pressure,zonal2.lat);
end
%Set up mesh grid
[X,Y,Z]=meshgrid(zonal2.lon,zonal2.lat,zonal2.depth);
%%
%Calculate Pressure Gradient between the two stations
%PG=(Pressure Difference)/Distance
%Calculate distance
[pdistance,az]=distance(zonal1.lat,zonal1.lon,zonal2.lat,zonal2.lon);
%calculate pressure difference
for i=1:length(zonal1.pressure)
for j=1:length(zonal2.pressure)
pdiff=[pdiff,zonal1.pressure(1)-zonal2.pressure(j)]';
end
end

Risposte (1)

Hassaan
Hassaan il 3 Gen 2024
% Initialize zonal1 and zonal2 structures
zonal1 = struct();
zonal2 = struct();
% --- Zonal 1 Data Preparation ---
% Change directory to Zonal 1 data
cd '/Users/shaymaalali/Desktop/Physical_Oceanography/cchdo_320620170820';
% Read Files
direc = dir('*.nc');
% Read data for zonal1
for i = 1:length(direc)
zonal1.lat = ncread(direc(i).name,'latitude');
zonal1.lon = ncread(direc(i).name,'longitude');
zonal1.pressure = ncread(direc(i).name,'pressure');
zonal1.temp = ncread(direc(i).name,'temperature');
zonal1.salinity = ncread(direc(i).name,'salinity');
end
% Calculate density for zonal1
for j = 1:length(zonal1.pressure)
[zonal1.density,alpha, beta] = gsw_rho_alpha_beta_CT_exact(zonal1.salinity, zonal1.temp, zonal1.pressure);
end
% Calculate depth from Pressure for zonal1
zonal1.depth = Depth(zonal1.pressure, zonal1.lat);
% --- Zonal 2 Data Preparation ---
% Change directory to Zonal 2 data
cd '/Users/shaymaalali/Desktop/Physical_Oceanography/cchdo_74EQ20151206';
% Unzip Files
unzip('1_a05_74EQ20151206_nc_ctd.zip');
unzip('0_a05_74EQ20151206_ct1.zip');
direc = dir('*.nc');
% Read data for zonal2
for i = 1:length(direc)
zonal2.lat = ncread(direc(i).name,'latitude');
zonal2.lon = ncread(direc(i).name,'longitude');
zonal2.pressure = ncread(direc(i).name,'pressure');
zonal2.temp = ncread(direc(i).name,'temperature');
zonal2.salinity = ncread(direc(i).name,'salinity');
end
% Calculate density for zonal2
for j = 1:length(zonal2.pressure)
[zonal2.density,alpha, beta] = gsw_rho_alpha_beta_CT_exact(zonal2.salinity, zonal2.temp, zonal2.pressure);
end
% Calculate depth from Pressure for zonal2
zonal2.depth = Depth(zonal2.pressure, zonal2.lat);
% --- Calculate Pressure Gradient between the two stations ---
% Calculate the great-circle distance between zonal1 and zonal2
[pdistance, az] = distance(zonal1.lat, zonal1.lon, zonal2.lat, zonal2.lon);
% Convert distance from degrees to meters
distance_in_km = deg2km(pdistance);
distance_in_meters = distance_in_km * 1000;
% Initialize the pressure difference array
pdiff = [];
% Calculate pressure difference for each pair
for i = 1:length(zonal1.pressure)
for j = 1:length(zonal2.pressure)
pdiff = [pdiff; zonal1.pressure(i) - zonal2.pressure(j)];
end
end
% Calculate the pressure gradient
pressure_gradient = pdiff ./ distance_in_meters;
% Display results
disp('Pressure Gradient (Pa/m): ');
disp(pressure_gradient);
Notes:
  • Ensure the distance, deg2km, and gsw_rho_alpha_beta_CT_exact functions are available in your MATLAB environment. They might require specific toolboxes or external libraries.
  • This code assumes that the Depth function is defined elsewhere in your scripts or path. Make sure it's available.
  • Ensure that the lengths of zonal1.pressure and zonal2.pressure match or adjust the code accordingly if they differ.
  • The code calculates the pressure gradient between every possible pair of points from zonal1 and zonal2. If you have a large number of points, this could become computationally intensive.
  • Before running the script, verify the paths and function availability in your MATLAB environment.
------------------------------------------------------------------------------------------------------------------------------------------------
If you find the solution helpful and it resolves your issue, it would be greatly appreciated if you could accept the answer. Also, leaving an upvote and a comment are also wonderful ways to provide feedback.
Professional Interests
  • Technical Services and Consulting
  • Embedded Systems
  • Electrical and Electronics Engineering

Categorie

Scopri di più su Programming in Help Center e File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by