Best projection for Getis Ord Gi* Statistic Analysis

I am about to calculate the Getis Ord Gi* statistic for some rasters whose extents are global.

• In choosing a projection, do I want to preserve distance (e.g. equidistant) to obtain accurate results?

• What do people recommend as a decent global projection for such an analysis?

Gi* doesn't care about distance, it cares about weights (which you could set to inverse distance if you want...) so you need to think about how to form weights between raster cells on the whole globe.

This is tricky, because at the poles, your cells are a different shape and area to those on the equator, and have a different adjacency relationship with neighbouring cells. Indeed, if you have, say, a 1000x1000 resolution raster over the whole globe then you have 1000 thin cells all touching at a point at the N and S pole, but at the equator each cell has four adjacent neighbours that share a long border. The Gi* weights, which are ad-hoc at the best of times, would be hard to justify.

By ad-hoc I mean that there's no method for choosing the "best" weights. For a problem with small irregular areas, typically you'd choose adjacency and use a 1/0 weight, or you'd use proportion of common boundary, or distance, or... well, anything.

Common boundary proportion might be a good option. At the equator the four neighbouring cells share 1/4 of their boundary in common. At the poles, cells adjacent E-W share nearer half their boundary, and cells N-S share much smaller fractions. Across the pole they share zero, but then they are very thin there and have closer relations with their E-W neighbours than the ones across the top.

But Gi* is a bit exploratory anyway so as long as you clearly state your assumptions as regards weights I don't think you'll get much kickback. But following up with something more statistically rigourous (model-based) would be advisable..

• Can you recommend a software package in which you would perform this? Commented Apr 23, 2015 at 16:59
• The function `localG` in the R package `spdep` can do this, and gives you full freedom in the weights matrix. Doc in here: cran.r-project.org/web/packages/spdep/spdep.pdf Commented Apr 23, 2015 at 17:14

Proposal: Conic Equidistant Projection

This preserves distance over all meridians and two parallels of choice (could be optimized for major areas of interest). Notice, most distance distortion is in the southern oceans out of the area of interest.

As an alternative, you could divide the world into 5 identical projections centered on each continent. This creates a natural way to break up the job, teasing out regional hotspot emergence and then gluing it back together for global visualization

A couple of contributing factors:

• all of the spatial-relationship approaches used in hotspot analysis (Inverse distance, Zone of Indifference, etc) derive their weights based on planar euclidian distance
• a reality of global analysis for land qualities is that spatial significance does not cross giant oceans, so you could really be talking about a 5-continent approach
• If you capture the points from existing rasters, you will preserve density, however the original projection introduces varying densification of data into the Getis-Ord Gi* algorithm. I’m not sure how this affects the outcomes, but it might be that the original projection has more influence

I'm currently struggling with similar issues, so feel free to wait for a more authoritative answer. However, I can't recommend an equidistant projection. It sounds great, but keep in mind that the distances are true only along specific lines. Measured along other directions, the distances will not be true. I doubt that you want to use a single projected raster, because it won't take into account the fact that the earth wraps around from one end of the raster to the other.

So that's it for the less speculative portion of the "answer". I would go with a conformal projection, although that's probably not ideal for your purposes, either. If it's possible, you might consider chopping the raster into overlapping strips 3 UTM zones wide, projecting them to the central zone, running the analysis on each strip, then clipping to the central zone. This would give you a bunch of locally calculated statistics derived from moderately accurate data. [This idea has not been evaluated by me or the experts.]

If you really want to be accurate, you should probably write your own software to calculate the distance weights from a 3d surface.

PS Hope someone else posts saying, "Actually, this existing scientific software already does what you want."