Calculating area volume from longitude, latitude and altitude

Hello! I'm trying to calculate the volume of airspace(3d polygon) from its longitude, latitude and altitude.

The polygon's vertices are given with longitude and latitude (ex. [128 129 129 128 128], [37 37 38 38 37])

and its altitude ranging from 1000~5000 feet.

How can I get the volume of this 3d polygon?

Answer 1

[darova]

• Convert spherical coordinates into cartesian using sph2cart
• Use alphaShape to build an object and calculate volume

[Walter Roberson]

airspaces have curved edges. You cannot convert the endpoints to cartesian and calculate volumes from that.

You could sub-divide into smaller values that individually had smaller curvature and use cartesian coordinates on those.

Or you could calculate in spherical coordinates.

[hye wook Kim]

Thanks for the kind help for both @darova and @Walter Roberson :)

As the sph2cart function gets [azimuth, elevation, r] as input data, how can I convert latitude, longitude and altitude into those input form?

+@Walter Roberson Is there any formula calculating the volume with spherical coordinates?

[darova]

@Walter Roberson, indeed, thanks

[darova]

Something is missing: how to close polygon since altitude (height/radius) is different at each point?

az0 = [128 129 129 128 128];        % longitude
el0 = [37 37 38 38 37];             % latitude
R0 = linspace(1000,5000,numel(az0)-1);% altitude
R0 = [R0 R0(1)];                    % assume that first and last point have the same altitude
t0 = 1:numel(az0);                  % parameter (number of points)
    % create 100 points from 5
t1 = linspace(t0(1),t0(end),100);   
az1 = interp1(t0,az0,t1);           
el1 = interp1(t0,el0,t1);           
R1 = interp1(t0,R0,t1);
[x1,y1,z1] = sph2cart(az1,el1,R1); % convert ot cartesian
surf([x1;x1*0],[y1;y1*0],[z1;z1*0],'facecolor','r')
line(x1,y1,z1)
axis equal
light

[darova]

• how can I convert latitude, longitude and altitude into those input form?

[hye wook Kim]

Sorry for short explanation. The polygon that I want to shape is cuboid that has same lat and long for vertices.

I'm wordering how to chanege the code above to calcuate the vomule of cuboid:)

[darova]

Here is the way to plot it. But it's not a cube. Any idea how to calculate volume?

[az,el] = meshgrid(128:0.2:129, 37:0.2:38);

[x0,y0,z0] = sph2cart(az,el,1000); % convert ot cartesian
[x1,y1,z1] = sph2cart(az,el,5000); % convert ot cartesian
plot3([x0(:) x1(:)]',[y0(:) y1(:)]',[z0(:) z1(:)]','b')
surface(x0,y0,z0,'facecolor','g')
surface(x1,y1,z1,'facecolor','g')
axis equal
light

[hye wook Kim]

I was tied up with other work recently. Let me find and try other approaches and share it in this post:)

[Walter Roberson]

You need to change the 1000 and 5000 to "radius of the Earth in that area, in feet, plus this".

https://rechneronline.de/earth-radius/

Latitude at 37, Earth radius at sea level: 3958.404 miles

Latitude at 38, Earth radius at sea level: 3958.18 miles

That is a decrease of

format long g
(3958.404 - 3958.18) * 5280

ans = 
          1182.72000000085

a bit over 1000 feet. So 1182.72 feet at 38 degrees (which is within the airspace) would be at sea level at 37 degrees (and so not in the airspace.)

So, you cannot treat the base of the airspace as flat: the curvature difference is equal to approximately 1/4 of the difference in altitudes between the 1000 and 5000 foot level.

That page gives a radius equation that it might be useful to integrate over.

latitude B, radius R, radius at equator r1, radius at pole r2

R = √ [ (r1² * cos(B))² + (r2² * sin(B))² ] / [ (r1 * cos(B))² + (r2 * sin(B))² ]

Integrate altitude + R over altitude and latitude in radians, and multiply by radians spanned by longitude, to get volume.

[hye wook Kim]

Thank you for your help! @Walter Roberson.

As you said, I'm trying to solve this problem using the triplet integral (integral3 in matlab).

and I got 3 questions as follows;

Q1. How to calculate the volume in mile unit? (ex. 1000 nm3)

Q2. Does R need to be calculated in the loop for every latitude?

Q3. how to construct formula for integral3 to apply your recommand; 'Integrate altitude + R over altitude and latitude in radians, and multiply by radians spanned by longitude'.

This is the code that I wrote and seems to show wrong answer...

%% Coordinates(Lat,Lon,Alt) of Target Airspace
Latitude = [128 129 129 128 128]; % Degree
Longitude = [37 37 38 38 37]; % Degree
Bottom_Altitude = distdim(1000,'feet','nauticalmiles'); % Convert FT to NM
Top_Altitude = distdim(5000,'feet','nauticalmiles'); % Convert FT to NM
%% Get Centroid coordinates of Target Airspace
polyin = polyshape(Latitude, Longitude);
[centroid_latitutde, centroid_longitude] = centroid(polyin);
%% Get Radius of the Earth in target airspace
r1 = 3963.191; % Earth's radius (nm) at sea level at the EQUATOR.
r2 = 3949.903; % Earth's radius (nm) at sea level at the POLES.
B = centroid_latitutde; % Latitude of target airspace's centroid.
R = sqrt(((r1^2*cos(B))^2+(r2^2*sin(B))^2)/((r1*cos(B))^2+(r2*sin(B))^2)); % Earth's radius (nm) at sea level from Target centroid
%% Get Spherical Coordinates from Lat,Lon,Alt
min_rho = Bottom_Altitude + R; max_rho = Top_Altitude + R; % Height (nm) from the Earth's center.
min_phi = deg2rad(min(Latitude)); max_phi = deg2rad(max(Latitude)); % Latitude
min_theta = deg2rad(min(Longitude)); max_theta = deg2rad(max(Longitude)); % Longitude
%% Apply Triplet Integral to get the volume
vol_fun = @(rho,phi,theta) rho.*phi.*theta;
vol_Boundary = integral3(vol_fun,min_rho,max_rho,min_phi,max_phi,min_theta,max_theta);