8
$\begingroup$

I'm a big fan of the new Mathematica 10 geographic capabilities and functions such as GeoDisk[], GeoCircle[] and others. One limitation of these functions, however, is the transformation of the actual GeoGraphics primitives into a list of coordinates. For GeoDisk[] for instance I can always do this manually via functions such as:

Clear[genperimeter]
genperimeter[citypos_, radius_] := 
  Reverse@Normal@
      LatitudeLongitude[GeoDestination[citypos, {radius, #}]] & /@ 
   Most[Range[0, 360, 30]];

which simply calculates points spaced in 30 degree increments around the geographic circle. It would be much more convenient however to extract the perimeter coordinates directly from a GeoDisk[] object. I've spend some time looking through the documentation but have not found a solution.

Many thanks ahead of time for any tips.

Mac

$\endgroup$

1 Answer 1

1
$\begingroup$

Here's a way:

Let's use a GeoDisk[] centered around the great city of Houston, TX:

houmap = GeoGraphics@
GeoDisk[CityData[{"Houston", "Texas", "UnitedStates"}, "Position"], Quantity[50, "Miles"]]

enter image description here

Now we can extract all the drawn polygons in the image:

Cases[houmap, Polygon[__], Infinity] // Length

3

To verify which one we want:

Graphics/@ Cases[houmap, Polygon[__], Infinity]

enter image description here

So we want the last Polygon:

Cases[houmap, Polygon[__], Infinity][[-1, 1]]

{{-95.3863,32.0591},{-95.3571,32.0586},<<177>>,{-95.4156,32.0586},{-95.3863,32.0591}}

These points, however, are x-y points, not the true latitude longitude of the Disk. To convert, we have write a function:

findcoordinates[geographic_, polylocation_] := Module[
{pts = Cases[geographic, Polygon[__], Infinity][[polylocation, 1]], proj, geomodel},
{proj, geomodel} = {GeoProjection, GeoModel} /. 
 AbsoluteOptions[geographic, {GeoProjection, GeoModel}];
(GeoPosition[GeoGridPosition[#1, proj, geomodel]] & /@ pts)[[All, 1]]
]

Now we can apply the function:

findcoordinates[houmap, -1]

{{30.5064,-95.3863},{30.5059,-95.3571},<<177>>,{30.5059,-95.4156},{30.5064,-95.3863}}

We can verify by using a GeoPath:

GeoPath[findcoordinates[houmap, -1]] // GeoGraphics

enter image description here

$\endgroup$
4
  • 1
    $\begingroup$ You could just ask for {GeoProjection, GeoModel} /. AbsoluteOptions[geographics, {GeoProjection, GeoModel}] directly instead of tearing apart the coordinate display function. $\endgroup$ Oct 7, 2014 at 14:35
  • 1
    $\begingroup$ @BrettChampion Man, that's way easier. Edit in progress... $\endgroup$
    – kale
    Oct 7, 2014 at 14:42
  • $\begingroup$ Brilliant - many thanks. I'm wondering for cases whether the first two geographic objects will always be those related to the background display ? If this is true the object extraction process can be automated and findcoordinates[] will only require one parameter. $\endgroup$
    – Mac
    Oct 8, 2014 at 8:24
  • $\begingroup$ @Mac On the cases I tested, it was always third, but no promises. $\endgroup$
    – kale
    Oct 8, 2014 at 14:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.