Converting lat/long to x/y coordinates, to generate random variates from a Bivariate Normal Distribution

Question

I am trying to convert latitude and longitude points to x/y coordinates, generate random variates from a bivariate normal distribution, and then convert those random x/y points back into latitude longitude coordinates. I am trying to get them to fall between 3 standard deviations of the center point that I pick, but I don't think that it is working. Here is my code so far.

    XYPos = GeoGridPosition[GeoPosition[{46.6572, -41.3552}], "Mercator"]
    mu1mu2 = Extract[XYPos, {1}]
    RandomXYs = 
      RandomVariate[BinormalDistribution[mu1mu2, {30/3, 30/3}, 0], 1000]
    RandomLatLongs = GeoPosition[GeoGridPosition[RandomXYs, "Mercator"]]

I am trying to get the 1000 points that I generate to fall within a circle 30 nautical miles in radius, clearly my units are off. Any help would be appreciated!

Answer 1

Thirty nautical miles is about 1/2 degree, so your standard deviation should be 1/6 degree. Let's test with a small sample.

RandomXYs = RandomVariate[BinormalDistribution[mu1mu2, {1/6, 1/6}, 0], 4];
GeoPosition[GeoGridPosition[RandomXYs, "Mercator"]]

GeoPosition[
    {{46.6226, -41.4243}, {46.6964, -41.3399}, {46.4473, -41.4126}, {46.8159, -41.4419}}]

Looks OK to me.

Answer 2

I am not sure what the aim is here. I post this hoping it may be useful. Using data from another question.

data = Reverse /@ {{47.5164, -36.5325}, {38.4084, -44.8774}, {37.829, \
-26.4011}, {43.415, -30.1504}, {40.9715, -50.8068}, {58.4615, \
-42.449}, {34.7541, -33.7088}, {55.2115, -32.6238}, {53.9798, \
-44.9295}, {46.6124, -39.4725}, {46.6667, -34.707}, {52.2855, \
-43.9319}, {50.7799, -38.798}, {50.277, -30.7066}, {50.6877, \
-46.4707}, {41.2522, -36.6759}, {61.4961, -44.8129}, {45.7457, \
-40.0388}, {38.0716, -35.5958}, {46.8175, -49.1262}, {49.4128, \
-44.9692}, {40.5598, -42.7279}, {55.1421, -34.8597}, {55.2034, \
-44.0431}, {40.7811, -48.9856}, {47.2722, -58.9129}, {38.3931, \
-27.3422}, {48.8027, -29.0668}, {45.8126, -42.8723}, {44.9651, \
-41.559}, {45.1573, -35.1754}, {49.7213, -54.5828}, {42.8633, \
-30.4128}, {47.5673, -31.3586}, {40.5429, -32.6634}, {50.9218, \
-44.1694}, {46.1213, -50.7524}, {54.7713, -59.3347}, {52.6424, \
-29.7903}, {53.2222, -44.2971}, {45.0114, -52.9965}, {54.8556, \
-31.8349}, {22.4978, -35.2292}, {49.4626, -64.454}, {45.5471, \
-57.4164}, {48.219, -59.1224}};
Needs["MultivariateStatistics`"]
Manipulate[
 GeoGraphics[{Cases[
    Graphics[{Red, Opacity[0.4], EllipsoidQuantile[data, p]}], 
    Line[x__] :> {Red, Opacity[0.2], Polygon[x]}, Infinity], {White, 
    Point[data]}, 
   Text[Style["+", Red, 20], MultivariateTrimmedMean[data, 0]]}, 
  GeoBackground -> "ReliefMap", 
  Frame -> True], {p, {0.5, 0.75, 0.95, 0.99}}]