Finding number of points inside a 2.5*2.5 lat-lon grid box using matlab?
Show older comments
Hello everyone,
I have a 1-D longitude and latitude data (attached excel file) (info about the data: irregular spacing and not in a increasing/ decreasing order).
I want to find a number of this irregular lat-lon points coming on/in each grid box each of size 2.5*2.5 of a regular lat-lon grid.
Therefore, I created a regular lat lon meshgrid using,
T=readtable('Question_data.xlsx'); lo=T.LON; la=T.LAT; % irregular data
res = 2.5; % grid box size
lon = -180:res:180;
lat = -90:res:90;
[lon,lat] = meshgrid(lon,lat);
At the end I will have a matrix of dimension 73*145 (i.e. in regula grid of resolution 2.5) containing how many points coming on/in each grid box.
I know something can be done with knnsearch but not sure how to execute it defining no of points and distance. I hope you have understood my requirement.
Looking forward to your help. Thanks in advance.
Answers (1)
T=readtable('Question_data.xlsx');
res = 2.5; % grid box size
% lon = -180:res:180;
% lat = -90:res:90;
lon = linspace(-180,180,round(360/res+1));
lat = linspace(-90,90,round(180/res+1));
% bin the data in 2 dimensions:
N = histcounts2(T.LAT,T.LON,lat,lon); % 72-by-144 (number of bins is one less than the number of edges)
disp(N);
disp(nnz(N)); % number of bins that have stuff
6 Comments
Subhodh Sharma
on 18 Jun 2022
Edited: Subhodh Sharma
on 18 Jun 2022
What makes the most sense to me is to use the bin centers for the contour, rather than the bin edges. (And that has the added benefit that the sizes work out, since N has an element for each bin, and every bin has one center).
T=readtable('Question_data.xlsx');
res = 2.5; % grid box size
lon = linspace(-180,180,round(360/res+1));
lat = linspace(-90,90,round(180/res+1));
% bin the data in 2 dimensions:
N = histcounts2(T.LAT,T.LON,lat,lon); % 72-by-144 (number of bins is one less than the number of edges)
% make the contour according to the centers of the bins, not the edges
lon_centers = (lon(1:end-1)+lon(2:end))/2;
lat_centers = (lat(1:end-1)+lat(2:end))/2;
contourf(lon_centers,lat_centers,N,'EdgeColor','none')
% or maybe a nice surface instead:
% surface(lon_centers,lat_centers,N,'EdgeColor','none','FaceColor','interp')
xlim([-180 180]);
hold on
C = load('coast.mat');
plot(C.long,C.lat,'k')
Subhodh Sharma
on 21 Jun 2022
Voss
on 21 Jun 2022
As written there
rr=nansum(nansum(N));
rr is the total number of points in all N.
Of course you mean
rr=nansum(nansum(ind));
However, nansum is not recommended; you can use the 'omitnan' option in sum instead to achieve the same result:
rr=sum(sum(ind,'omitnan'),'omitnan');
And you can avoid having to call sum twice, either by converting ind into a vector:
rr=sum(ind(:),'omitnan');
or by using the 'all' option, if you have R2018b or later:
rr=sum(ind,'all','omitnan');
Of course, in this case, since you are summing elements of a matrix that comes from histcounts2, I wouldn't expect there to be any NaNs at all (though I haven't done any tests to see under what conditions histcounts2 returns a NaN count, if any), so I would expect that you can probably safely omit 'omitnan', e.g.:
rr=sum(ind(:));
Finally, ind is defined with the lat and lon idx swapped. It should be:
ind=N(idx_lat,idx_lon);
(Notice N has T.LAT first in histcounts2.)
N = histcounts2(T.LAT,T.LON,lat,lon);
Subhodh Sharma
on 21 Jun 2022
Edited: Subhodh Sharma
on 21 Jun 2022
@Voss thanks again!
I don't know whether this below statement of yours is right or wrong?
Because if I try to calculate ind using
ind=N(idx_lat,idx_lon);
It throws an error "Index exceeds matrix dimensions" and this is obviously expected.
I guess the calculated N first dimension represents lon and second dimension lat. I also checked with reversed condition,
N = histcounts2(T.LON,T.LAT,lon,lat);
here N has dimension of 144 by 72.
lon has dimension of 145 and lat has dimension of 73. If the statement "number of bins is one less than the number of edges" is true. then the logic has to consistent.
Correct me if I am wrong. Thanks in advance.
Let me try it both ways.
All this stuff doesn't depend on how N is defined, so I do it once at the top:
T=readtable('Question_data.xlsx');
res = 2.5; % grid box size
lon = linspace(-180,180,round(360/res+1));
lat = linspace(-90,90,round(180/res+1));
lat_centers = (lat(1:end-1)+lat(2:end))/2;
lon_centers = (lon(1:end-1)+lon(2:end))/2;
idx_lat = lat_centers > 0 & lat_centers < 40;
idx_lon = lon_centers > 110 & lon_centers < 150;
First case: defining N with lat first, lon second:
% lat first, lon second:
N = histcounts2(T.LAT,T.LON,lat,lon); % 72-by-144
size(N)
% lat first, lon second:
ind=N(idx_lat,idx_lon) % no error here
Second case: defining N with lon first, lat second:
% lon first, lat second:
N = histcounts2(T.LON,T.LAT,lon,lat); % 144-by-72
size(N)
% lon first, lat second:
ind=N(idx_lon,idx_lat) % no error here
So I'm unable to reproduce the error "Index exceeds matrix dimensions". You must be doing something different than one of these two ways.
To address your second point:
"here N has dimension of 144 by 72.
lon has dimension of 145 and lat has dimension of 73. If the statement "number of bins is one less than the number of edges" is true. then the logic has to consistent."
- lon and lat are bin edges.
- lon_centers and lat_centers are bin centers.
- N is counts inside the bins.
Therefore, N has size [one less than the length of lat]-by-[one less than the length of lon] (or the reverse, depending on which way you define N).
Another way to say the same: N has size [length of lat_centers]-by-[length of lon_centers] (or the reverse, depending on which way you define N).
Another way to say the same: In a given dimension (lat or lon), each bin has two edges and one center. N is bin counts, so there's one element of N for each bin, i.e., one count for each bin center.
Categories
Find more on JSON Format in Help Center and File Exchange
on 16 Jun 2022
on 21 Jun 2022
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
