How to get longitude and latitude?

Question

Ara 2024-5-27

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2122876-how-to-get-longitude-and-latitude

评论： Ara 2024-5-29

Dear All,

I want to calculate longtitude and latitude for total electron content (TEC) using GNSS radio occulttaion (RO) data. I have the X-GPS, Y-GPS, and Z-GPS as well as X-LEO, Y-LEO, and Z-LEO. How can I calculate it in Matlab?

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Manikanta Aditya 2024-5-27

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2122876-how-to-get-longitude-and-latitude#answer_1463701

在 MATLAB Online 中打开

Hi @Ara

To calculate the longitude and latitude from the given GNSS (Global Navigation Satellite System) radio occultation data, you can transform the coordinates from the Earth-Centered, Earth-Fixed (ECEF) coordinate system to the Geodetic coordinate system (latitude, longitude, and altitude).

Here’s a simple example of how you might do this:

function [lat, lon, alt] = ecef2lla_custom(x, y, z)
    % Constants for the WGS-84 ellipsoid
    a = 6378137.0; % semi-major axis in meters
    f = 1/298.257223563; % flattening
    e2 = f * (2-f); % square of eccentricity
    % Calculations
    lon = atan2(y, x);
    p = sqrt(x.^2 + y.^2);
    theta = atan2(z * a, p * (1 - f) * a);
    lat = atan2(z + e2 * (1 - f)^2 * sin(theta).^3 * a, p - e2 * cos(theta).^3 * a);
    N = a ./ sqrt(1 - e2 * sin(lat).^2);
    alt = p ./ cos(lat) - N;
    % Convert to degrees
    lat = rad2deg(lat);
    lon = rad2deg(lon);
end
% Example coordinates
X_GPS = 1123456.7;
Y_GPS = 2123456.7;
Z_GPS = 3123456.7;
X_LEO = 2123456.7;
Y_LEO = 3123456.7;
Z_LEO = 4123456.7;
% Call the custom function to convert GNSS ECEF coordinates to geodetic coordinates
[lat_GPS, lon_GPS, alt_GPS] = ecef2lla_custom(X_GPS, Y_GPS, Z_GPS);
% Call the custom function to convert LEO ECEF coordinates to geodetic coordinates
[lat_LEO, lon_LEO, alt_LEO] = ecef2lla_custom(X_LEO, Y_LEO, Z_LEO);
% Display the results for GPS
fprintf('GPS Coordinates:\n');
GPS Coordinates:
fprintf('Latitude: %.6f degrees\n', lat_GPS);
Latitude: 52.733047 degrees
fprintf('Longitude: %.6f degrees\n', lon_GPS);
Longitude: 62.118037 degrees
fprintf('Altitude: %.2f meters\n', alt_GPS);
Altitude: -2424368.20 meters
% Display the results for LEO
fprintf('\nLEO Coordinates:\n');
LEO Coordinates:
fprintf('Latitude: %.6f degrees\n', lat_LEO);
Latitude: 47.728614 degrees
fprintf('Longitude: %.6f degrees\n', lon_LEO);
Longitude: 55.790493 degrees
fprintf('Altitude: %.2f meters\n', alt_LEO);
Altitude: -774830.00 meters

I hope this helps!

32 个评论
显示 30更早的评论隐藏 30更早的评论

Manikanta Aditya 2024-5-27

编辑：Manikanta Aditya 2024-5-27

在 MATLAB Online 中打开

@Ara

Follow this:

Download and read the real data from https://data.cosmic.ucar.edu/gnss-ro/cosmic2/provisional/spaceWeather/level2/2023/001/
Process the electron density profile to calculate Slant TEC (STEC).
Convert the Slant TEC (STEC) to Vertical TEC (VTEC).

% Download and load the NetCDF file
filename = 'ionPrf_C002.2023001.00.00.G26_20232331922.27.nc'; 
filepath = fullfile('path_to_downloaded_file', filename); 
% Read the data from the NetCDF file
altitudes = ncread(filepath, 'MSL_alt'); % Altitude in km
electronDensity = ncread(filepath, 'elec_dens'); % Electron density in electrons per cubic meter
% Ensure data is sorted by altitude in descending order
[altitudes, idx] = sort(altitudes, 'descend');
electronDensity = electronDensity(idx);
% Calculate Slant TEC (STEC) by integrating electron density
STEC = trapz(altitudes, electronDensity) * 1e-5; % Convert to TEC units (TECU, 1 TECU = 10^16 electrons/m^2)
% Function to convert STEC to VTEC
function VTEC = convert_stec_to_vtec(STEC, elevationAngle)
    % Mapping function for Single Layer Model (SLM)
    mappingFunction = 1 / sqrt(1 - (cosd(elevationAngle)^2 * (6371 / (6371 + 350))^2)); % Example with ionospheric height of 350 km
    VTEC = STEC / mappingFunction;
end
% Example elevation angle (replace with actual angle if available)
elevationAngle = 30; % Example elevation angle in degrees
% Convert STEC to VTEC
VTEC = convert_stec_to_vtec(STEC, elevationAngle);
% Display the results
fprintf('Total Electron Content (STEC): %.2f TECU\n', STEC);
fprintf('Vertical Total Electron Content (VTEC): %.2f TECU\n', VTEC);

Hope this helps.

Manikanta Aditya 2024-5-27

在 MATLAB Online 中打开

I think, this should answer your requirement query:

% Define the file path and file name (replace with actual file path and name)
filename = 'ionPrf_C002.2023001.00.00.G26_20232331922.27.nc'; 
filepath = fullfile('path_to_downloaded_file', filename); 
% Read the data from the NetCDF file
altitudes = ncread(filepath, 'MSL_alt'); % Altitude in km
electronDensity = ncread(filepath, 'elec_dens'); % Electron density in electrons per cubic meter
latitudes = ncread(filepath, 'lat'); % Latitude in degrees
longitudes = ncread(filepath, 'lon'); % Longitude in degrees
% Ensure data is sorted by altitude in descending order
[altitudes, idx] = sort(altitudes, 'descend');
electronDensity = electronDensity(idx);
% Calculate Slant TEC (STEC) by integrating electron density
STEC = trapz(altitudes, electronDensity) * 1e-5; % Convert to TEC units (TECU, 1 TECU = 10^16 electrons/m^2)
% Function to convert STEC to VTEC
function VTEC = convert_stec_to_vtec(STEC, elevationAngle)
    % Mapping function for Single Layer Model (SLM)
    mappingFunction = 1 / sqrt(1 - (cosd(elevationAngle)^2 * (6371 / (6371 + 350))^2)); % Example with ionospheric height of 350 km
    VTEC = STEC / mappingFunction;
end
% Example elevation angle 
elevationAngle = 30; 
% Convert STEC to VTEC
VTEC = convert_stec_to_vtec(STEC, elevationAngle);
% Display the results
fprintf('Total Electron Content (STEC): %.2f TECU\n', STEC);
fprintf('Vertical Total Electron Content (VTEC): %.2f TECU\n', VTEC);
% Plot VTEC as a function of geographic longitude and latitude
figure;
scatter(longitudes, latitudes, 100, VTEC, 'filled');
colorbar;
title('VTEC as a Function of Geographic Longitude and Latitude');
xlabel('Longitude (degrees)');
ylabel('Latitude (degrees)');
c = colorbar;
c.Label.String = 'VTEC (TECU)';

Ara 2024-5-27

在 MATLAB Online 中打开

Thank You. Please do not upset with me. I understand that you are providing your great help. I simply do not know how to do it so I ask my questions. I simp-ly need your help. You wrote a great code so I wish to get the plot of it.

I untar all the files in this folder (C:\Users\I7-2600K\Documents\Paper\Thesis). that is my file path. Is that correct?

I got file name from one of the file that is this "ionPrf_C2E1.2023.001.00.00.R19_0001"

filename = 'ionPrf_C2E1.2023.001.00.00.R19_0001'; 
filepath = fullfile('C:\Users\I7-2600K\Documents\Paper\Thesis', filename);
%also correct this part
altitudes = ncread(filepath, 'MSL_alt'); % Altitude in km
electronDensity = ncread(filepath, 'ELEC_dens'); % Electron density in electrons per cubic meter
latitudes = ncread(filepath, 'lat'); % Latitude in degrees
longitudes = ncread(filepath, 'lon'); % Longitude in degrees
But I receive same error.
Error using internal.matlab.imagesci.nc/openToRead (line 1299)
Unable to open file "C:\Users\I7-2600K\Documents\Paper\Thesis\ionPrf_C2E1.2023.001.00.00.R19_0001".
Error in internal.matlab.imagesci.nc (line 124)
                    this.openToRead();
Error in ncread (line 71)
ncObj   = internal.matlab.imagesci.nc(ncFile);
Error in MainTest2 (line 6)
altitudes = ncread(filepath, 'MSL_alt'); % Altitude in km

Manikanta Aditya 2024-5-27

在 MATLAB Online 中打开

The error message indicates that the input color argument for the scatter function is incorrect. This can happen if VTEC contains invalid values or if it's not in the correct format.

To address this issue, we need to ensure that VTEC contains valid values and that the color input for scatter is properly formatted.

% Plot VTEC as a function of geographic longitude and latitude
figure;
scatter(longitudes, latitudes, 100, VTEC, 'filled');
colorbar;
title('VTEC as a Function of Geographic Longitude and Latitude');
xlabel('Longitude (degrees)');
ylabel('Latitude (degrees)');
c = colorbar;
c.Label.String = 'VTEC (TECU)';

To create a contour plot of VTEC, you can use the contour or contourf functions. Here's an example of how to create a contour plot using contourf:

% Create a meshgrid for longitude and latitude
[X, Y] = meshgrid(longitudes, latitudes);
% Reshape VTEC to match the dimensions of the meshgrid
VTEC_matrix = reshape(VTEC, size(X));
% Create the contour plot
figure;
contourf(X, Y, VTEC_matrix);
colorbar;
title('Contour Plot of VTEC');
xlabel('Longitude (degrees)');
ylabel('Latitude (degrees)');

Manikanta Aditya 2024-5-27

在 MATLAB Online 中打开

Negative values for TEC are not physically meaningful and indicate that there might be an issue with the data processing.

filename = 'ionPrf_C2E1.2023.001.00.00.R19_0001.nc'; 
filepath = fullfile('C:\Users\I7-2600K\Documents\Paper\Thesis', filename);
% Verify that the file exists
if ~isfile(filepath)
    error('File not found: %s', filepath);
end
% Read the data from the NetCDF file
altitudes = ncread(filepath, 'MSL_alt'); % Altitude in km
electronDensity = ncread(filepath, 'ELEC_dens'); % Electron density in electrons per cubic meter
latitudes = ncread(filepath, 'lat'); % Latitude in degrees
longitudes = ncread(filepath, 'lon'); % Longitude in degrees
% Check for any negative electron density values and correct them if necessary
if any(electronDensity < 0)
    warning('Negative electron density values found. Setting them to zero.');
    electronDensity(electronDensity < 0) = 0;
end
% Ensure data is sorted by altitude in ascending order for integration
[altitudes, idx] = sort(altitudes);
electronDensity = electronDensity(idx);
% Calculate Slant TEC (STEC) by integrating electron density
STEC = trapz(altitudes, electronDensity) * 1e-5; % Convert to TEC units (TECU, 1 TECU = 10^16 electrons/m^2)
% Function to convert STEC to VTEC
function VTEC = convert_stec_to_vtec(STEC, elevationAngle)
    % Mapping function for Single Layer Model (SLM)
    mappingFunction = 1 / sqrt(1 - (cosd(elevationAngle)^2 * (6371 / (6371 + 350))^2)); % Example with ionospheric height of 350 km
    VTEC = STEC / mappingFunction;
end
% Example elevation angle (replace with actual angle if available)
elevationAngle = 30; % Example elevation angle in degrees
% Convert STEC to VTEC
VTEC = convert_stec_to_vtec(STEC, elevationAngle);
% Verify VTEC values
if any(isnan(VTEC)) || any(isinf(VTEC))
    error('VTEC contains NaN or Inf values. Please check your data.');
end
% Ensure longitudes, latitudes, and VTEC have the same length
if length(longitudes) ~= length(latitudes) || length(latitudes) ~= length(VTEC)
    error('longitudes, latitudes, and VTEC must have the same length.');
end
% Plot VTEC as a function of geographic longitude and latitude
figure;
scatter(longitudes, latitudes, 100, VTEC, 'filled');
colorbar;
title('VTEC as a Function of Geographic Longitude and Latitude');
xlabel('Longitude (degrees)');
ylabel('Latitude (degrees)');
c = colorbar;
c.Label.String = 'VTEC (TECU)';
% Create a meshgrid for longitude and latitude
[lonGrid, latGrid] = meshgrid(unique(longitudes), unique(latitudes));
% Interpolate VTEC onto the grid
VTEC_grid = griddata(longitudes, latitudes, VTEC, lonGrid, latGrid, 'natural');
% Create the contour plot
figure;
contourf(lonGrid, latGrid, VTEC_grid);
colorbar;
title('Contour Plot of VTEC');
xlabel('Longitude (degrees)');
ylabel('Latitude (degrees)');

Manikanta Aditya 2024-5-27

在 MATLAB Online 中打开

Yes, it is possible to set negative values to zero and ensure that the lengths of longitudes, latitudes, and VTEC are the same. You can do this by using interpolation or simply filtering out the data points that do not match.

filename = 'ionPrf_C2E1.2023.001.00.00.R19_0001.nc'; % Add the file extension ".nc"
filepath = fullfile('C:\Users\I7-2600K\Documents\Paper\Thesis', filename);
% Verify that the file exists
if ~isfile(filepath)
    error('File not found: %s', filepath);
end
% Read the data from the NetCDF file
altitudes = ncread(filepath, 'MSL_alt'); % Altitude in km
electronDensity = ncread(filepath, 'ELEC_dens'); % Electron density in electrons per cubic meter
latitudes = ncread(filepath, 'lat'); % Latitude in degrees
longitudes = ncread(filepath, 'lon'); % Longitude in degrees
% Check for any negative electron density values and correct them if necessary
electronDensity(electronDensity < 0) = 0;
% Ensure data is sorted by altitude in ascending order for integration
[altitudes, idx] = sort(altitudes);
electronDensity = electronDensity(idx);
% Calculate Slant TEC (STEC) by integrating electron density
STEC = trapz(altitudes, electronDensity) * 1e-5; % Convert to TEC units (TECU, 1 TECU = 10^16 electrons/m^2)
% Function to convert STEC to VTEC
function VTEC = convert_stec_to_vtec(STEC, elevationAngle)
    % Mapping function for Single Layer Model (SLM)
    mappingFunction = 1 / sqrt(1 - (cosd(elevationAngle)^2 * (6371 / (6371 + 350))^2)); % Example with ionospheric height of 350 km
    VTEC = STEC / mappingFunction;
end
% Example elevation angle (replace with actual angle if available)
elevationAngle = 30; % Example elevation angle in degrees
% Convert STEC to VTEC
VTEC = convert_stec_to_vtec(STEC, elevationAngle);
% Ensure all three arrays have the same length
minLength = min([length(longitudes), length(latitudes), length(VTEC)]);
longitudes = longitudes(1:minLength);
latitudes = latitudes(1:minLength);
VTEC = VTEC(1:minLength);
% Verify that all arrays have the same length after filtering
if length(longitudes) ~= length(latitudes) || length(latitudes) ~= length(VTEC)
    error('longitudes, latitudes, and VTEC must have the same length.');
end
% Plot VTEC as a function of geographic longitude and latitude
figure;
scatter(longitudes, latitudes, 100, VTEC, 'filled');
colorbar;
title('VTEC as a Function of Geographic Longitude and Latitude');
xlabel('Longitude (degrees)');
ylabel('Latitude (degrees)');
c = colorbar;
c.Label.String = 'VTEC (TECU)';
% Create a meshgrid for longitude and latitude
[lonGrid, latGrid] = meshgrid(unique(longitudes), unique(latitudes));
% Interpolate VTEC onto the grid
VTEC_grid = griddata(longitudes, latitudes, VTEC, lonGrid, latGrid, 'natural');
% Create the contour plot
figure;
contourf(lonGrid, latGrid, VTEC_grid);
colorbar;
title('Contour Plot of VTEC');
xlabel('Longitude (degrees)');
ylabel('Latitude (degrees)');

Ara 2024-5-28

在 MATLAB Online 中打开

Thank you very much for your great explanations.

I asked ChatGPT and it mentioned "Typical DCB values for GNSS RO satellites and receivers are in the range of 1 to 20 nanoseconds, depending on the specific frequencies and systems involved. " I got 10nanosecond for satellite and 5 nanosecond for receivers. And then change it to TECU and it gives me this range. Does it sound correct to you? In one day we should have more points not only a point. Do you have any idea to get the desirable figures? The contour plot is not working? It mentioned X should be scalar but how we can do it?

satelliteDCB = 61.73; % DCB for the satellite in TECU
receiverDCB = 30.86;

Ara 2024-5-29

Thank you, Manikanta Aditya.

You helped me alot.

All the best,

Ara

请先登录，再进行评论。

How to get longitude and latitude?

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

32 个评论
显示 30更早的评论隐藏 30更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

How to get longitude and latitude?

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

32 个评论 显示 30更早的评论隐藏 30更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

WeChat

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

32 个评论
显示 30更早的评论隐藏 30更早的评论