3
$\begingroup$
GeoGraphics[{
    {
        GeoStyling["StreetMap"],
        First[GeoNearest["Country",geoPosition]]["Polygon"],GeoDisk[geoPosition,Quantity[300,"Miles"]]
    }
        GeoMarker[geoPosition]
},GeoBackground->None]

Suppose I evaluate the above expression with

geoPosition = GeoPosition[{35.6762`,139.6503`}]

It plots the background with a union of country's polygon and disk's polygon but I want intersection of the two. How to do this?

$\endgroup$
3
  • $\begingroup$ I have more-or-less done this in this answer. I simply used RegionIntersection after making Regions from the polygons. $\endgroup$
    – Carl Lange
    Feb 22, 2020 at 19:42
  • $\begingroup$ @CarlLange do you know how to get polygon structure for the GeoDisk[geoPosition,Quantity[300,"Miles"]]? $\endgroup$
    – user13892
    Feb 22, 2020 at 20:11
  • $\begingroup$ Easiest thing to do is generate a number of points equidistant from the centre and create a polygon from that. I'll add an answer in a moment. $\endgroup$
    – Carl Lange
    Feb 22, 2020 at 20:40

1 Answer 1

6
$\begingroup$

We can do this fairly easily by converting our geographic regions into Regions and using RegionIntersection.

First let's get our geographic regions.

centre = GeoPosition[{35.6762`, 139.6503`}]

p1 = First[GeoNearest["Country", centre]]["Polygon"]

Now, it's not easy to use GeoDisk directly for our calculation, so we'll regenerate this by creating a Polygon with a list of points equidistant from the centre.

radius = Quantity[300, "Miles"]

p2 = Polygon[
  Table[GeoDestination[centre, GeoDisplacement[{radius, b}]], {b, 0, 360, 
    1}]]

Now we will double-check that this looks right:

GeoGraphics[{p1, p2}]

enter image description here

Now we compute the intersection by turning the polygons into regions and taking RegionIntersection of that.

r1 = Region[p1 /. GeoPosition[x___] -> x]

r2 = Region[p2 /. GeoPosition[x___] -> x]

int = RegionIntersection[r1, r2, PerformanceGoal -> "Speed"]

enter image description here

Now we can turn this back into a geo polygon:

geoint = int /. Region[Polygon[x_]] -> Polygon[GeoPosition[x]]

And there we are:

GeoGraphics[{geoint}]

enter image description here

Now we can do things like GeoArea easily with GeoArea[geoint] (in this case, 178393km^2!)

Note that Sato island is not included - you may need to use GeoVariant[..., "AllAreas"] to get every extended part of a country. It may return FilledCurves which may be harder to turn into Regions.

$\endgroup$
3
  • $\begingroup$ One more question, doesn't RegionIntersection also return a polygon object if the inputs are polygon? What is the purpose of MeshPrimitives[BoundaryDiscretizeRegion[...], 2]? $\endgroup$
    – user13892
    Feb 22, 2020 at 20:58
  • $\begingroup$ Indeed, thank you, I have edited my answer. That was a holdover from some previous code I used as the base of this answer. $\endgroup$
    – Carl Lange
    Feb 22, 2020 at 21:03
  • $\begingroup$ Sorry I am wrong RegionIntersection returns MeshRegion in the case of polygon inputs so I think MeshPrimitives[..., 2] is necessary to get 2D-polygons. $\endgroup$
    – user13892
    Feb 22, 2020 at 21:59

Your Answer

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

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