Plot for Month and year

Dear All,

I have data with format of NETCDF file. I can read one of them and plot it. For each day I have many data so I wish to read all of them to provide result for one day and then expand it to provide it for month and year. Also, I would like to have a data in .mat format, for example for each month. Would you please help me in this regard?

Best,

Ara

[Manikanta Aditya]

Hi @Ara, if possible can you share the NETCDF file here!

[Ara]

Hi Manikanta Aditya,

I could not share the data as it mentioned the format is not accepted. Here is the link for data to provide a better view of it.

https://data.cosmic.ucar.edu/gnss-ro/cosmic2/provisional/spaceWeather/level2/2023/

[Manikanta Aditya]

Okay thanks!

Answer 1

[Manikanta Aditya]

Hi, @Ara,

Based on the understanding from your explaination in the question, you can use some MATLAB functions and script for it.

Use can use the 'ncread' function: https://in.mathworks.com/help/matlab/ref/ncread.html

Check it out:

% Define paths to your NETCDF files
netcdfPath = 'path/to/netcdf/files/';

% Define output directory for processed data
outputDir = 'path/to/output/directory/';

% Loop through each NETCDF file
netcdfFiles = dir(fullfile(netcdfPath, '*.nc'));
for fileIdx = 1:numel(netcdfFiles)
    % Read NETCDF file
    ncFile = fullfile(netcdfPath, netcdfFiles(fileIdx).name);
    % Read necessary data from NETCDF file (use appropriate commands for your data)
    data = ncread(ncFile, 'variable_name');
    % Process daily data (e.g., calculate daily averages)
    dailyAverage = mean(data, 'all');
    
    dailyDataFileName = fullfile(outputDir, ['daily_data_', datestr(netcdfFiles(fileIdx).datenum, 'yyyymmdd'), '.mat']);
    save(dailyDataFileName, 'dailyAverage');
end

% Aggregate daily data to monthly
monthlyDataFileName = fullfile(outputDir, ['monthly_data_', datestr(netcdfFiles(fileIdx).datenum, 'yyyymm'), '.mat']);
save(monthlyDataFileName, 'monthlyData');

% Aggregate monthly data to yearly
yearlyDataFileName = fullfile(outputDir, ['yearly_data_', datestr(netcdfFiles(fileIdx).datenum, 'yyyy'), '.mat']);
save(yearlyDataFileName, 'yearlyData');

• Consider this workaround script as a base and it should help you to start off.

I hope this gives you some help with your query!

[Ara]

Thank you for the code.

I have received an error as below.

Error using datestr (line 187)
Unable to convert input to specified date string.
Failed to convert input to date numbers.
Error in MainPath (line 22)
monthlyDataFileName = fullfile(outputDir, ['monthly_data_', datestr(netcdfFiles(fileIdx).datenum, 'yyyymm'), '.mat']);

[Manikanta Aditya]

Interesting! I feel might be an issue with extracting the date from the NETCDF file names.

% Loop through each NETCDF file
netcdfFiles = dir(fullfile(netcdfPath, '*.nc'));
for fileIdx = 1:numel(netcdfFiles)
    % Extract date from the file name
    fileName = netcdfFiles(fileIdx).name;
    % Assuming the date is in the format YYYYMMDD, extract it from the file name
    fileDate = fileName(1:8); % Adjust this line according to the actual format of your file names
    
    % Read NETCDF file
    ncFile = fullfile(netcdfPath, fileName);
    % Read necessary data from NETCDF file 
    data = ncread(ncFile, 'variable_name');
    
    % Process daily data (e.g., calculate daily averages)
    dailyAverage = mean(data, 'all');
    
    % Save processed daily data
    dailyDataFileName = fullfile(outputDir, ['daily_data_', fileDate, '.mat']);
    save(dailyDataFileName, 'dailyAverage');
    
    monthYear = fileDate(1:6); % Extract YYYYMM from the file name
    monthlyDataFileName = fullfile(outputDir, ['monthly_data_', monthYear, '.mat']);

    year = fileDate(1:4); % Extract YYYY from the file name
    yearlyDataFileName = fullfile(outputDir, ['yearly_data_', year, '.mat']);
end

In this code, we extract the date part from the file name (fileDate) assuming that it's in the format YYYYMMDD. Adjust the extraction logic according to the actual format of your file names. Then, we use this date to construct the file names for saving the processed data.

Please try to fix the issue with your end as well! @Ara

[Ara]

Perfect.It works perfectly for daily. THANK You!

[Manikanta Aditya]

Great to know @Ara!

[Ara]

I want to read the daily mat file and then contour plot VTEC. How can I match it with this code? I can load the mat file and then...?...

% 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, 'GEO_lat'); % Latitude in degrees
longitudes = ncread(filepath, 'GEO_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
% Example DCB values (replace with actual values)
satelliteDCB = 61.73; % DCB for the satellite in TECU
receiverDCB = 30.86; 
% Adjust STEC for DCB
STEC_adjusted = STEC - (satelliteDCB + receiverDCB);
% Convert adjusted STEC to VTEC
VTEC_adjusted = convert_stec_to_vtec(STEC_adjusted, elevationAngle);
% 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]

@Ara,

To create a contour plot of Vertical Total Electron Content (VTEC) from the daily .mat files, you'll need to follow these steps:

1. Load Data from .mat Files
2. Interpolate and Contour Plot

% Define the path to your .mat files
matPath = 'path/to/output/directory/';
matFiles = dir(fullfile(matPath, 'daily_data_*.mat'));

% Initialize arrays to store the combined data
allLongitudes = [];
allLatitudes = [];
allVTEC = [];

% Loop through each .mat file and load the data
for fileIdx = 1:numel(matFiles)
    % Load the daily data file
    matFile = fullfile(matPath, matFiles(fileIdx).name);
    data = load(matFile);
    
    % Assuming the data structure contains longitudes, latitudes, and VTEC
    longitudes = data.longitudes;
    latitudes = data.latitudes;
    VTEC = data.VTEC;
    
    % Append the data to the combined arrays
    allLongitudes = [allLongitudes; longitudes];
    allLatitudes = [allLatitudes; latitudes];
    allVTEC = [allVTEC; VTEC];
end

% Ensure all three arrays have the same length
minLength = min([length(allLongitudes), length(allLatitudes), length(allVTEC)]);
allLongitudes = allLongitudes(1:minLength);
allLatitudes = allLatitudes(1:minLength);
allVTEC = allVTEC(1:minLength);

% Create a meshgrid for longitude and latitude
[lonGrid, latGrid] = meshgrid(unique(allLongitudes), unique(allLatitudes));

% Interpolate VTEC onto the grid
VTEC_grid = griddata(allLongitudes, allLatitudes, allVTEC, 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)');
c = colorbar;
c.Label.String = 'VTEC (TECU)';

[Ara]

Thank you for the code and explanation.

I have to calculate VTEC from STEC. I have electron density that I need to read from NETCDF file and then convert it to VTEC. And then plot it as a contourplot.