Risk of generating sampling points from Pseudo-Mercator projection

Question

Is there any risk in using the Pseudo Mercator projection (EPSG:3857) to generate lat/long coordinates for locating sampling plots with a handheld GPS?

My understanding is that since it uses WGS-84 datum, there should be no issue -- but am also unsure as EPSG:3857 is an unrecognized projection and there might be idiosyncrasies I'm unaware of.

Is it necessary to reproject to EPSG:4326 or a UTM based SRS before deriving sampling point lat/long coordinates?

Answer 1

"No issue?" Hardly. Any time you generate random data you should do a chi-squared test to verify the randomness of your generator. Your distribution will be skewed toward the north or south poles, because the values were uniformly generated in Mercator (which infinitely exaggerates areas at the poles).

If you use an equal area projection with a uniform distribution, your points will be uniformly distributed by area. If you use a geographic coordinate system, you need to apply a trig function to the latitude value to counter the smaller area as the latitude approaches the poles.

This Python code demonstrates my point:

import math
from random import uniform

ubins = [0] * 15
tbins = [0] * 15
for i in range(0,1000):
    uval = uniform(0.0,90.0)
    ubins[int(math.floor(uval/6.0))] += 1
    fval = uval / 90.0
    tval = math.degrees(math.asin(fval))
    tbins[int(math.floor(tval/6.0))] += 1

print "\n   bin\trange\tuniform\ttrig"
for i in range(0,15):
    print '   {:d}\t{:2.0f}-{:2.0f}\t{:d}\t{:d}'.format(
                i,6.0*i,6.0*(i+1),ubins[i],tbins[i])

And the generated output:

C:\Temp>python distrib.py

   bin  range   uniform trig
   0     0- 6   64      91
   1     6-12   51      111
   2    12-18   82      113
   3    18-24   65      85
   4    24-30   71      110
   5    30-36   60      98
   6    36-42   72      83
   7    42-48   83      72
   8    48-54   65      66
   9    54-60   74      44
   10   60-66   68      44
   11   66-72   63      40
   12   72-78   55      22
   13   78-84   68      15
   14   84-90   59      6

The first column is the 6 degree bin number, with the range specified in the second column. The third column is a uniform distribution of values in the bin (by latitude), and the fourth column shows a uniform distribution by area on a sphere.

This, of course, only distributes values across the northern hemisphere, but the same procedure can be applied to the globe:

C:\Temp>python globe.py

   bin  range   uniform trig
   0    -90,-78 69      13
   1    -78,-66 71      26
   2    -66,-54 64      50
   3    -54,-42 72      72
   4    -42,-30 67      81
   5    -30,-18 68      96
   6    -18,-6  62      113
   7    -6, 6   63      115
   8     6,18   76      84
   9    18,30   72      92
   10   30,42   59      79
   11   42,54   64      72
   12   54,66   63      63
   13   66,78   57      37
   14   78,90   73      7

The code to generate this is left as an exercise, but I will warn that had to use arc-cosine instead of arc-sine, because the starting bin was at the pole, not the equator.

All in all, I'd recommend a simple uniform distribution on an appropriate equal-area projection, then deproject to get geographic WGS84 values.