Skip to content

Density analysis producing extreme values?

An answer to this question on the GIS Stack Exchange.

Question

I have point data of gas flares and I am trying to create a simple density raster of those flares. Results are either 0 or in the 100,000-1,000,000 range (when there's 5 to 40 or so points in each box).

At least as I see it, I am trying to do the most basic of point density analyses. I'm just simply trying to count up how many points are in that box! I'm using a reference raster of 15 arcseconds. My point density dialogue is as below. Note that I've tried it two different ways with very similar results. Any suggestions of where I'm going wrong? I'm guessing it has to do with my search radius, but how do I tell Arc I want my search radius to be the same as the pixel?

Dialogue: tool dialog

Results: results

Environments: Processing extent same as reference raster
Input: Gas Flares Points
Population: none
Output: ...\Flare_Dens.tif
Output cell size: 0.004167 (15 arcsecs in dec degrees, same as raster I'm snapping to)
Neighborhood: Rectangle
Radius: Map, 0.004167 x 0.004167

OR (getting similar results in the >100,000 range)
Radius: Cell, 1 x 1

Answer

The problem with the projections others have suggested is that you'll be unable to divide the spherical surface of the Earth into equally-sized cells using a rectangular grid. A discrete global grid would be most appropriate here as it will generate a grid of equally-sized, hexagonal cells.

I don't think Arc has the capability to do such things, but my package dggridR makes it easy enough if you have access to R:

library(dggridR)
library(dplyr)
#Construct a global grid with cells approximately 1000 miles across
dggs          <- dgconstruct(spacing=1000, metric=FALSE, resround='down')
#Load included test data set
data(dgquakes)
#Get the corresponding equally-sized, hexagonal grid cells for each 
#earthquake epicenter (lat-long pair)
dgquakes$cell <- dgtransform(dggs,dgquakes$lat,dgquakes$lon)
#Get the number of earthquakes in each equally-sized cell
quakecounts   <- dgquakes %>% group_by(cell) %>% summarise(count=n())

Plotting the result yields:

Binning on a discrete global grid using dggridR, as project onto a Mercator projection

The cells do not appear to be of equal size because of the distortion imposed by the projection. The method others have suggested is to do such a projection and lay a rectalinear grid on top of it. From the image above, it should be obvious what a bad idea that can be. Choosing projections carefully and working in limited geographical areas reduces this error, but only so much.

If we reproject the same data onto a sphere, it is obvious that the cells are of equal size, as desired:

Binning on a discrete global grid using dggridR, as projected onto a sphere

Naturally, you could choose a finer grid as appropriate for your data set.