How to make a lightning strokes/ square km / day map

I have WWLLN data in .loc format. I have 92 files, from where I need to select the lats and lons (strokes) for +/-30 location only. How can I select this range from 92 .loc files using a loop and how to make a density map ( Lightning Strokes / square km/ day) using Matlab R2017a?

Answer 1

[Gautam]

Hello Joydeb,

You can follow the below workflow to read files, select the required data points and create the density map

1. Loop through each file, read the data, and extract the latitude and longitude values.
2. Select only the data points that fall within the specified latitude and longitude range (±30 degrees)
3. Store the filtered data for all files
4. Create density map using the “imagesc” function

Here’s a sample code that demonstrates this workflow

% Directory containing .loc files
dataDir = 'path_to_your_loc_files';
files = dir(fullfile(dataDir, '*.loc'));

% Initialize arrays to store latitude and longitude data
allLats = [];
allLons = [];

% Loop through each .loc file
for k = 1:length(files)
    % Construct full file path
    filePath = fullfile(dataDir, files(k).name);
    
    % Read the .loc file
fid = fopen(filePath, 'r');
    data = textscan(fid, '%f%f%*[^\n]', 'Delimiter', ' ', 'CollectOutput', true);
    fclose(fid);
    
    % Extract latitude and longitude
    latLonData = data{1}; % Assuming latitudes and longitudes are in the first two columns
    lats = latLonData(:, 1);
    lons = latLonData(:, 2);
    
    % Filter data for ±30 range
    validIdx = abs(lats) <= 30 & abs(lons) <= 30;
    allLats = [allLats; lats(validIdx)];
    allLons = [allLons; lons(validIdx)];
end

% Define the grid
lat_min = min(allLats);
lat_max = max(allLats);
lon_min = min(allLons);
lon_max = max(allLons);

% Define the resolution of the grid (e.g., 0.1 degree)
lat_res = 0.1;
lon_res = 0.1;

% Create the grid
lat_edges = lat_min:lat_res:lat_max;
lon_edges = lon_min:lon_res:lon_max;

% Count the number of strokes in each grid cell
counts = histcounts2(latitudes, longitudes, lat_edges, lon_edges);

% Calculate the area of each grid cell in square kilometers
% Approximation: 1 degree latitude ~ 111 km, 1 degree longitude ~ 111 km * cos(latitude)
lat_km = 111; % km per degree latitude
lon_km = @(lat) 111 * cosd(lat); % km per degree longitude

% Calculate the area of each grid cell
[lat_grid, lon_grid] = meshgrid(lat_edges(1:end-1) + lat_res/2, lon_edges(1:end-1) + lon_res/2);
cell_area_km2 = lat_km * lat_res * lon_km(lat_grid) * lon_res;
% Calculate density (strokes per square km per day)
% Assume you have data for a certain number of days, e.g., 30 days
num_days = 30;
density = counts ./ (cell_area_km2 * num_days);

% Plot the density map
figure;
imagesc(lon_edges, lat_edges, density);
set(gca, 'YDir', 'normal'); % Correct the Y-axis direction
colorbar;
xlabel('Longitude');
ylabel('Latitude');

You can refer to the following documentations for more information on the MATLAB functions used:

1. textscan: https://www.mathworks.com/help/matlab/ref/textscan.html
2. imagesc; https://www.mathworks.com/help/matlab/ref/imagesc.html