# Method to generate points in any projection

I am writing an application in ArcMap to generate a grid of points. These points have a very defined layout and I have worked out an algorithm that works great under a projected coordinate system. However, as soon as I try it under a geographic coordinate system, my math doesn't stand up as well and I am not an expert in this type of math (but will to learn if pointed in the right direction)

Are there any suggested methods, libraries, or anything else that might make this easier?

A little bit more information about the dots themselves: In a circle of 144 foot radius, the dots will be generated in a hex type pattern. The hex is made up of an equilateral triangle with sides of 9 meters (in meters because ArcMap IPoint X,Y properties are in meters).

• Can you get the inverse projection? It seems like you could use that to preserve your distribution. Oct 28, 2010 at 22:43
• Wait are you asking (in a different form) for a projection from spherical coordinates to a flat surface and back that preserves length and area? Because I'm pretty sure the answer is no. Oct 28, 2010 at 22:49
• @Jared #2 In a way I guess I am asking for that. But would it be possible to do it in a projected coordinate system then project it or would that also break down? Oct 28, 2010 at 23:03
• @Jared There's no requirement, not even implicitly, for area to be preserved. (That doesn't matter anyway: no projection even preserves length.) Moreover, with a remarkably small number of exceptions (corresponding to the Platonic solids) there's no such thing as a regular grid on the earth's surface anyway: so we should be thinking in terms of approximations from the outset. Oct 29, 2010 at 14:07

@om_henners' suggestion will work well, especially if you are creating the grid in the "right" way.

The right way views each grid point as being offset from an origin by the sum of integral multiples of two basis vectors. In projected coordinates a good choice of basis vectors for a triangular grid consists of one arbitrary vector of the desired length (e.g., 9 meters) together with a rotation of that vector by 120 degrees. For example, if you want the bases of the triangles to be horizontal you can pick the first vector to be e = (9, 0). The second one would then equal f = (-9/2, Sqrt(3)*9/2) = (-4.5, 7.794229). Letting O = (Ox, Oy) be the coordinates of the origin, all other points on the grid are of the form

i * e + j * f + O

where i and j are integers. Here's the punchline: to compensate for the distortion in geographic coordinates, divide the first coordinate of each basis vector by the cosine of the latitude, convert all offsets to degrees, but otherwise everything else remains the same.

For example, at latitude 40 degrees the cosine is 0.766044. This converts e to e' = (9/0.766044, 0) = (11.74867, 0) and f becomes f' = (-4.5/0.766044, 7.794229) = (-5.874333, 7.794229). Now convert these four coordinates to degrees: division by 111,111 (meters per degree) is usually close enough. Using these vectors instead of e and f in the grid formula, but keeping O at the same location (and using the same sets of i and j values) lays out your grid in geographic coordinates. Over these short distances (up to 144 feet in this case, or out to several kilometers in general) you can get extremely high accuracy by first adjusting the latitude to a spherically-equivalent latitude before computing the cosine, but typically the error is so small (much less than 1%) that it doesn't matter.

The vector approach is illustrated in an old article of mine.

BTW, a nice feature is that you can easily rotate and shift the grid. To rotate it, just rotate the original basis vectors e and f by the desired amount, adjust their first coordinates by the cosine, and then recompute the grid points as linear combinations of these adjusted rotated basis vectors. To shift the grid, just shift the origin O. Thus you are no longer a slave to the coordinate system: you can orient and position the grid as needed for the application. Another nice feature is that every mathematically possible regular grid can be constructed this way. Just change e and f. For example, when f is at 90 degrees to e you will create a rectangular grid.

@baens Yeah, I agree with @Jared - There's no way with a projection to a sphere and back again to preserve length and area. You could potentially use one of the "compromise" projections (see wikipeda), but otherwise you may have to choose.

That said, at your scale you could use the quick and dirty 111,111 metres is one degree in the y direction, and 111,111 metres * cos (latitude) is one degree in the x direction. This is from @whuber 's answer to offsetting a lat/long by metres.

Is there a reason why you wouldn't create the features in a SRS that is true for distance and area and then project/transform the features to a geographic SRS when you are done?

Maybe something like Matt Perry's Tissot Indicatrix example? http://www.perrygeo.net/wordpress/?p=4 Features are created in an orthographic projection and then projected to geographic.

I guess the answer really depends on what you are trying to accomplish. If you want the points to reflect the same point on the ground, you will want to create them in a different SRS and project them to geographic. If you don't really care to preserve the location of each point, but want the shape to appear true in geographic, I would create them directly there.

• Those Tissot Indicatrices are not quite right. They will be ok in the center where distortion is low but (in general) they cannot even be computed near the limits of the projection's domain. The approach is a good one but it needs modification. What works well in practice is to generate true geodetic circles of about 10 - 100 meters in radius, project those, and then uniformly rescale them around their centers (in map coordinates) to make them visible. Additionally, a deeper analysis of this procedure lets you compute useful projection properties: see quantdec.com/tissot . Oct 29, 2010 at 14:04