The following code works, but is inefficient. I would like to find a better implementation.
The problem is the following. I have a list of points, representing a matrix of 2 columns and N rows. Each point represents an event occurring on a map. I create a tessellation for the given map, creating a grid of x-by-y squares I want to calculate the number of events that fall into each square.
The result should be an x-by-y matrix where an element gives the number of events falling inside a given square. I will use this to do some more point pattern analysis, etc.
grid
is created as a list of polygons, that correspond to the squares of the tessellation. It is created as follows, with x0
and y0
defining the “origin corner” of the map from which to start the tessellation:
grid =
Table[
Polygon[
{{x0 + (i - 1)*dsq, y0 + (j - 1)*dsq},
{x0 + i*dsq, y0 + (j - 1)*dsq},
{x0 + i*dsq, y0 + j*dsq},
{x0 + (i - 1)*dsq, y0 + j*dsq},
{x0 + (i - 1)*dsq, y0 + (j - 1)*dsq}}],
{i, 1, xGrid}, {j, 1, yGrid}];
{{Polygon[{{41.645, -87.884}, {41.6828, -87.884}, {41.6828, -87.8463}, {41.645, -87.8463}, {41.645, -87.884}}], Polygon[{{41.645, -87.8463}, {41.6828, -87.8463}, {41.6828, -87.8085}, {41.645, -87.8085}, {41.645, -87.8463}}], Polygon[{{41.645, -87.8085}, {41.6828, -87.8085}, {41.6828, -87.7708}, {41.645, -87.7708}, {41.645, -87.8085}}], ...}}
I initialize the matrix to store the results (resMat
)
xGrid = 10;
yGrid = 10;
resMat = Table[0, {x, 1, xGrid}, {y, 1, yGrid}] ;
This is the expression I use to go through each square and calculate how many points are in it:
Table[
If[RegionMember[grid[[x, y]], points[[n, All]]],
resMat[[x, y]] = resMat [[x, y]] + 1],
{n, 1, Length[points]}, {x, 1, xGrid}, {y, 1, yGrid}] ;
I use the following to quickly visualize the results:
ArrayPlot[Reverse[resMat], ColorFunction -> ColorData["SolarColors"]]
The computation of resMat
works, but is really slow and inefficient. The RegionMember[…]
is the bottleneck: for a grid size of 100-by-100, the analysis of a single point can take 10-20 sec. My point list contains hundreds of points, making the computation with the current code too time consuming.
I could not come up with a better implementation so far. Therefore, I will be glad to get some help.
Nearest
applied to the midpoints of the grid cells with optionDistanceFunction -> ChessboardDistance
to find all events within each grid cell (distance <= than half the cells' edge length. This should be much faster. $\endgroup$ – Henrik Schumacher Jun 1 '18 at 12:48