The first step is to convert lat,lon to some projected x,y values because lats and lons are not equally spaced. If you have the Mapping Toolbox you can use projfwd for the coordinate transformation.
When you have the coordinates in x,y space, then it's just a matter of creating a 2D histogram to get lightning strike density. You can search the File Exchange site for 2D histogram functions that others have already created. Or you could manually create a 2D histogram by creating a grid, then in a loop you could determine the sum of x,y values that are within the bounds of each grid cell. Here's an example:
x = 20*randn(300,1);
y = 10*randn(300,1)+x/3;
plot(x,y,'rx')
axis equal;
gridres = 10;
xgrid = -100:gridres:100;
ygrid = -50:gridres:50;
% Preallocate a histogram:
hist2 = NaN(length(ygrid),length(xgrid));
for row = 1:length(ygrid)-1
for col = 1:length(xgrid)-1
% Get indices of all x,y points inside this grid cell:
ind = x>=xgrid(col)-gridres/2 & x<xgrid(col+1)-gridres/2 & y>=ygrid(row)-gridres/2 & y<ygrid(row+1)-gridres/2;
%
% Log the sum of number of lighting strikes inside the grid cell:
hist2(row,col) = sum(ind);
end
end
%
hold on
h = imagesc(xgrid,ygrid,hist2/(gridres^2));
uistack(h,'bottom')
